diff --git a/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp b/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp index e3a13f6..34322d9 100644 --- a/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp +++ b/Modules/DiffusionCmdApps/Quantification/DiffusionIndices.cpp @@ -1,193 +1,196 @@ /*=================================================================== 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 "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #include #include #include /** * */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Diffusion Indices"); parser.setCategory("Diffusion Related Measures"); parser.setDescription("Computes requested diffusion related measures"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::String, "Input:", "input image (tensor, ODF or SH-coefficient image)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("index", "idx", mitkDiffusionCommandLineParser::String, "Index:", "index (fa, gfa, ra, ad, rd, ca, l2, l3, md, adc)", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string index = us::any_cast(parsedArgs["index"]); std::string outFileName = us::any_cast(parsedArgs["o"]); std::string ext = itksys::SystemTools::GetFilenameLastExtension(outFileName); if (ext.empty()) outFileName += ".nii.gz"; try { // load input image mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Diffusion Weighted Images", "SH Image", "ODF Image", "Tensor Image"}, {}); auto input = mitk::IOUtil::Load(inFileName, &functor); bool is_odf = (dynamic_cast(input.GetPointer()) || dynamic_cast(input.GetPointer())); bool is_dt = dynamic_cast(input.GetPointer()); bool is_dw = mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(input); if (is_odf) MITK_INFO << "Input is ODF image"; else if (is_dt) MITK_INFO << "Input is tensor image"; else if (is_dw) MITK_INFO << "Input is dMRI"; else { MITK_WARN << "Input is no ODF, SH, tensor or raw dMRI."; return EXIT_FAILURE; } mitk::LocaleSwitch localeSwitch("C"); if( is_odf && index=="gfa" ) { typedef itk::Vector OdfVectorType; typedef itk::Image OdfVectorImgType; OdfVectorImgType::Pointer itkvol; if (dynamic_cast(input.GetPointer())) + { + MITK_INFO << "Assuming MITK/MRtrix style SH convention!"; itkvol = mitk::convert::GetItkOdfFromShImage(input); + } else itkvol = mitk::convert::GetItkOdfFromOdfImage(input); typedef itk::DiffusionOdfGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(itkvol); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(gfaFilter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if( is_dt ) { typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(input.GetPointer()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itk_dti.GetPointer() ); if(index=="fa") measurementsCalculator->SetMeasure(MeasurementsType::FA); else if(index=="ra") measurementsCalculator->SetMeasure(MeasurementsType::RA); else if(index=="ad") measurementsCalculator->SetMeasure(MeasurementsType::AD); else if(index=="rd") measurementsCalculator->SetMeasure(MeasurementsType::RD); else if(index=="ca") measurementsCalculator->SetMeasure(MeasurementsType::CA); else if(index=="l2") measurementsCalculator->SetMeasure(MeasurementsType::L2); else if(index=="l3") measurementsCalculator->SetMeasure(MeasurementsType::L3); else if(index=="md") measurementsCalculator->SetMeasure(MeasurementsType::MD); else { MITK_WARN << "No valid diffusion index for input image (tensor image) defined"; return EXIT_FAILURE; } measurementsCalculator->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(measurementsCalculator->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else if(is_dw && (index=="adc" || index=="md")) { typedef itk::AdcImageFilter< short, double > FilterType; auto itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(input); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkVectorImagePointer ); filter->SetGradientDirections( mitk::DiffusionPropertyHelper::GetGradientContainer(input) ); filter->SetB_value( static_cast(mitk::DiffusionPropertyHelper::GetReferenceBValue(input)) ); if (index=="adc") filter->SetFitSignal(true); else filter->SetFitSignal(false); filter->Update(); itk::ImageFileWriter< itk::Image >::Pointer fileWriter = itk::ImageFileWriter< itk::Image >::New(); fileWriter->SetInput(filter->GetOutput()); fileWriter->SetFileName(outFileName); fileWriter->Update(); } else std::cout << "Diffusion index " << index << " not supported for supplied file type."; } catch (const itk::ExceptionObject& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp b/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp index 11c8ba6..02436a7 100644 --- a/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp +++ b/Modules/DiffusionCmdApps/Tractography/GlobalTractography.cpp @@ -1,133 +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()); 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 996fc2c..e6daff8 100644 --- a/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionCmdApps/Tractography/StreamlineTractography.cpp @@ -1,566 +1,569 @@ /*=================================================================== 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; RF)", us::Any(), false); + 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; 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); 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)); 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)); 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 == "Tensor" && params->m_Mode == mitk::StreamlineTractographyParameters::MODE::PROBABILISTIC)) + 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)); 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]; if (dynamic_cast(input.GetPointer())) { MITK_INFO << "Converting SH to ODF image"; - mitk::Image::Pointer mitkImg = dynamic_cast(input.GetPointer()); + 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); 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/Algorithms/itkShToOdfImageFilter.cpp b/Modules/DiffusionCore/Algorithms/itkShToOdfImageFilter.cpp index 1243cd1..2786812 100644 --- a/Modules/DiffusionCore/Algorithms/itkShToOdfImageFilter.cpp +++ b/Modules/DiffusionCore/Algorithms/itkShToOdfImageFilter.cpp @@ -1,80 +1,80 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkShToOdfImageFilter_cpp #define __itkShToOdfImageFilter_cpp #include #include #include #include "itkShToOdfImageFilter.h" #include #include namespace itk { template< class PixelType, int ShOrder > ShToOdfImageFilter< PixelType, ShOrder >::ShToOdfImageFilter() - : m_Toolkit(Toolkit::MRTRIX) + : m_Toolkit(mitk::ShImage::SH_CONVENTION::MRTRIX) { } template< class PixelType, int ShOrder > void ShToOdfImageFilter< PixelType, ShOrder >::BeforeThreadedGenerateData() { CalcShBasis(); } template< class PixelType, int ShOrder > void ShToOdfImageFilter< PixelType, ShOrder >::ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType) { typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); typename InputImageType::Pointer inputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); typedef ImageRegionConstIterator< InputImageType > InputIteratorType; InputIteratorType it(inputImage, outputRegionForThread); typedef ImageRegionIterator< OutputImageType > OutputIteratorType; OutputIteratorType oit(outputImage, outputRegionForThread); while(!it.IsAtEnd()) { auto pix = it.Get(); vnl_vector coeffs = pix.GetVnlVector(); OutputPixelType odf = (m_ShBasis * coeffs).data_block(); oit.Set(odf); ++it; ++oit; } } // generate spherical harmonic values of the desired order for each input direction template< class PixelType, int ShOrder > void ShToOdfImageFilter< PixelType, ShOrder >::CalcShBasis() { vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); - if (m_Toolkit==Toolkit::MRTRIX) + if (m_Toolkit==mitk::ShImage::SH_CONVENTION::MRTRIX) m_ShBasis = mitk::sh::CalcShBasisForDirections(ShOrder, U->as_matrix()); else m_ShBasis = mitk::sh::CalcShBasisForDirections(ShOrder, U->as_matrix(), false); } } #endif // __itkShToOdfImageFilter_cpp diff --git a/Modules/DiffusionCore/Algorithms/itkShToOdfImageFilter.h b/Modules/DiffusionCore/Algorithms/itkShToOdfImageFilter.h index dea8d8a..13e2dca 100644 --- a/Modules/DiffusionCore/Algorithms/itkShToOdfImageFilter.h +++ b/Modules/DiffusionCore/Algorithms/itkShToOdfImageFilter.h @@ -1,87 +1,83 @@ /*========================================================================= Program: Insight Segmentation & Registration Toolkit Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $ Language: C++ Date: $Date: 2006-03-27 17:01:06 $ Version: $Revision: 1.12 $ Copyright (c) Insight Software Consortium. All rights reserved. See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef __itkShToOdfImageFilter_h_ #define __itkShToOdfImageFilter_h_ #include #include +#include namespace itk{ /** \class ShToOdfImageFilter */ template< class PixelType, int ShOrder > class ShToOdfImageFilter : public ImageToImageFilter< itk::Image< itk::Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder >, 3 >, itk::Image< itk::Vector,3> > { public: - enum Toolkit { ///< SH coefficient convention (depends on toolkit) - FSL, - MRTRIX - }; - typedef itk::Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder > InputPixelType; typedef itk::Image InputImageType; typedef typename InputImageType::RegionType InputImageRegionType; typedef itk::Vector OutputPixelType; typedef itk::Image OutputImageType; typedef typename OutputImageType::RegionType OutputImageRegionType; typedef ShToOdfImageFilter Self; typedef itk::ImageToImageFilter Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(ShToOdfImageFilter, ImageToImageFilter) - itkSetMacro( Toolkit, Toolkit) ///< define SH coefficient convention (depends on toolkit) - itkGetMacro( Toolkit, Toolkit) ///< SH coefficient convention (depends on toolkit) + itkSetMacro( Toolkit, mitk::ShImage::SH_CONVENTION) ///< define SH coefficient convention (depends on toolkit) + itkGetMacro( Toolkit, mitk::ShImage::SH_CONVENTION) ///< SH coefficient convention (depends on toolkit) void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); protected: ShToOdfImageFilter(); ~ShToOdfImageFilter(){} void CalcShBasis(); vnl_matrix m_ShBasis; - Toolkit m_Toolkit; + mitk::ShImage::SH_CONVENTION m_Toolkit; private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkShToOdfImageFilter.cpp" #endif #endif //__itkShToOdfImageFilter_h_ diff --git a/Modules/DiffusionCore/Algorithms/itkShToRgbImageFilter.h b/Modules/DiffusionCore/Algorithms/itkShToRgbImageFilter.h index 771a112..2eb6dc9 100644 --- a/Modules/DiffusionCore/Algorithms/itkShToRgbImageFilter.h +++ b/Modules/DiffusionCore/Algorithms/itkShToRgbImageFilter.h @@ -1,154 +1,160 @@ /*=================================================================== 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 __itkShToRgbImageFilter_h #define __itkShToRgbImageFilter_h #include "itkUnaryFunctorImageFilter.h" #include "itkOrientationDistributionFunction.h" #include "itkRGBAPixel.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionIterator.h" #include +#include namespace itk { #define __IMG_DAT_ITEM__CEIL_ZERO_ONE__(val) (val) = \ ( (val) < 0 ) ? ( 0 ) : ( ( (val)>1 ) ? ( 1 ) : ( (val) ) ); /** \class ShToRgbImageFilter * */ template ,3> > class ShToRgbImageFilter : public ImageToImageFilter { public: /** Standard class typedefs. */ typedef ShToRgbImageFilter Self; typedef ImageToImageFilter Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef typename TInputImage::PixelType InputPixelType; typedef typename InputPixelType::ValueType InputValueType; /** Run-time type information (and related methods). */ itkTypeMacro( ShToRgbImageFilter, ImageToImageFilter ) /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Print internal ivars */ void PrintSelf(std::ostream& os, Indent indent) const override { this->Superclass::PrintSelf( os, indent ); } #ifdef ITK_USE_CONCEPT_CHECKING /** Begin concept checking */ itkConceptMacro(InputHasNumericTraitsCheck, (Concept::HasNumericTraits)); /** End concept checking */ #endif + itkSetMacro( ShConvention, mitk::ShImage::SH_CONVENTION) ///< define SH coefficient convention + itkGetMacro( ShConvention, mitk::ShImage::SH_CONVENTION) ///< SH coefficient convention + protected: ShToRgbImageFilter(){} ~ShToRgbImageFilter() override{} + mitk::ShImage::SH_CONVENTION m_ShConvention = mitk::ShImage::SH_CONVENTION::MRTRIX; + void ThreadedGenerateData( const typename OutputImageType::RegionType &outputRegionForThread, ThreadIdType) override { typename InputImageType::Pointer coeff_image = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); outputImage->SetSpacing( coeff_image->GetSpacing() ); // Set the image spacing outputImage->SetOrigin( coeff_image->GetOrigin() ); // Set the image origin outputImage->SetDirection( coeff_image->GetDirection() ); // Set the image direction outputImage->SetRegions( coeff_image->GetLargestPossibleRegion()); outputImage->Allocate(); typedef ImageRegionConstIterator< InputImageType > OdfImageIteratorType; OdfImageIteratorType it(coeff_image, outputRegionForThread); typedef ImageRegionIterator< OutputImageType > OutputImageIteratorType; OutputImageIteratorType oit(outputImage, outputRegionForThread); it.GoToBegin(); oit.GoToBegin(); typedef itk::OrientationDistributionFunction OdfType; - vnl_matrix sh2Basis = mitk::sh::CalcShBasisForDirections(ShOrder, itk::PointShell >::DistributePointShell()->as_matrix()); + vnl_matrix sh2Basis = mitk::sh::CalcShBasisForDirections(ShOrder, itk::PointShell >::DistributePointShell()->as_matrix(), m_ShConvention); while(!it.IsAtEnd() && !oit.IsAtEnd()) { InputPixelType x = it.Get(); Vector< float, ODF_SAMPLING_SIZE > odf_vals; vnl_vector< float > coeffs(x.GetNumberOfComponents()); for(unsigned int i=0; i dir; int pd = odf.GetPrincipalDiffusionDirectionIndex(); if (pd==-1) dir.fill(0); else dir = OdfType::GetDirection(pd); const float fa = odf.GetGeneralizedFractionalAnisotropy(); float r = fabs(dir[0]) * fa; float g = fabs(dir[1]) * fa; float b = fabs(dir[2]) * fa; float a = fa; __IMG_DAT_ITEM__CEIL_ZERO_ONE__(r); __IMG_DAT_ITEM__CEIL_ZERO_ONE__(g); __IMG_DAT_ITEM__CEIL_ZERO_ONE__(b); __IMG_DAT_ITEM__CEIL_ZERO_ONE__(a); OutputPixelType out; out.SetRed( r * 255.0f); out.SetGreen( g * 255.0f); out.SetBlue( b * 255.0f); out.SetAlpha( a * 255.0f); oit.Set(out); ++it; ++oit; } } private: ShToRgbImageFilter(const Self&); //purposely not implemented void operator=(const Self&); //purposely not implemented }; } // end namespace itk #endif diff --git a/Modules/DiffusionCore/IODataStructures/mitkShImage.cpp b/Modules/DiffusionCore/IODataStructures/mitkShImage.cpp index 90d273e..4fcf0ac 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkShImage.cpp +++ b/Modules/DiffusionCore/IODataStructures/mitkShImage.cpp @@ -1,144 +1,156 @@ /*=================================================================== 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 "mitkShImage.h" #include "mitkImageCast.h" #include "itkImage.h" #include "mitkImageVtkAccessor.h" #include #include "itkShToRgbImageFilter.h" mitk::ShImage::ShImage() : Image() , m_ShOrder(0) , m_NumCoefficients(0) + , m_ShConvention(SH_CONVENTION::MRTRIX) { m_RgbImage = nullptr; // not needed anymore as soon as all diffusion images are identified via properties anyway this->SetProperty("IsShImage", mitk::BoolProperty::New(true)); } mitk::ShImage::~ShImage() { } unsigned int mitk::ShImage::NumCoefficients() { if (m_NumCoefficients==0 || m_ShOrder==0) { m_NumCoefficients = static_cast(this->m_ImageDescriptor->GetChannelTypeById(0).GetNumberOfComponents()); int c=3, d=2-2*static_cast(m_NumCoefficients); int D = c*c-4*d; if (D>0) { int s = (-c+static_cast(sqrt(D)))/2; if (s<0) s = (-c-static_cast(sqrt(D)))/2; m_ShOrder = static_cast(s); } else if (D==0) m_ShOrder = static_cast(-c/2); } return m_NumCoefficients; } unsigned int mitk::ShImage::ShOrder() { NumCoefficients(); return m_ShOrder; } vtkImageData* mitk::ShImage::GetVtkImageData(int t, int n) { if(m_RgbImage.IsNull()) { ShOrder(); NumCoefficients(); ConstructRgbImage(); } return m_RgbImage->GetVtkImageData(t,n); } const vtkImageData*mitk::ShImage::GetVtkImageData(int t, int n) const { if(m_RgbImage.IsNull()) { ConstructRgbImage(); } return m_RgbImage->GetVtkImageData(t,n); } template void mitk::ShImage::Construct() const { typedef itk::Image,3> ImageType; typedef itk::ShToRgbImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); typename ImageType::Pointer itkvol = ImageType::New(); mitk::CastToItkImage(this, itkvol); filter->SetInput(itkvol); + filter->SetShConvention(this->m_ShConvention); filter->Update(); itk::Image,3>::Pointer tmp = filter->GetOutput(); m_RgbImage = mitk::Image::New(); m_RgbImage->InitializeByItk( tmp.GetPointer() ); m_RgbImage->SetVolume( tmp->GetBufferPointer() ); } void mitk::ShImage::ConstructRgbImage() const { switch (this->m_ImageDescriptor->GetChannelTypeById(0).GetNumberOfComponents()) { case 6: Construct<2>(); break; case 15: Construct<4>(); break; case 28: Construct<6>(); break; case 45: Construct<8>(); break; case 66: Construct<10>(); break; case 91: Construct<12>(); break; default : mitkThrow() << "SH order larger 12 not supported"; } } const vtkImageData* mitk::ShImage::GetNonRgbVtkImageData(int t, int n) const { return Superclass::GetVtkImageData(t,n); } vtkImageData* mitk::ShImage::GetNonRgbVtkImageData(int t, int n) { return Superclass::GetVtkImageData(t,n); } void mitk::ShImage::PrintSelf(std::ostream &os, itk::Indent indent) const { os << indent << "Spherical harmonics order: " << m_ShOrder << std::endl; os << indent << "Number of coefficients: " << m_NumCoefficients << std::endl; Superclass::PrintSelf(os, indent); } + +mitk::ShImage::SH_CONVENTION mitk::ShImage::ShImage::GetShConvention() const +{ + return m_ShConvention; +} + +void mitk::ShImage::ShImage::SetShConvention(SH_CONVENTION ShConvention) +{ + m_ShConvention = ShConvention; +} diff --git a/Modules/DiffusionCore/IODataStructures/mitkShImage.h b/Modules/DiffusionCore/IODataStructures/mitkShImage.h index 9aed144..34f9d61 100644 --- a/Modules/DiffusionCore/IODataStructures/mitkShImage.h +++ b/Modules/DiffusionCore/IODataStructures/mitkShImage.h @@ -1,71 +1,80 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkShImage__h #define __mitkShImage__h #include "mitkImage.h" #include "itkVectorImage.h" #include "mitkImageVtkAccessor.h" #include namespace mitk { /** * \brief this class encapsulates spherical harmonics coefficients */ class MITKDIFFUSIONCORE_EXPORT ShImage : public Image { public: + enum SH_CONVENTION { ///< SH coefficient convention (depends on toolkit) + FSL, // used in FSL, Dipy + MRTRIX // used in MRtrix, MITK + }; + typedef itk::Image ShOnDiskType; // we store the sh images in MRtrix 4D float format and convert on load to 3D multi-component images mitkClassMacro( ShImage, Image ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) virtual const vtkImageData* GetNonRgbVtkImageData(int t = 0, int n = 0) const; virtual vtkImageData* GetNonRgbVtkImageData(int t = 0, int n = 0); const vtkImageData* GetVtkImageData(int t = 0, int n = 0) const override; vtkImageData* GetVtkImageData(int t = 0, int n = 0) override; virtual void ConstructRgbImage() const; unsigned int ShOrder(); unsigned int NumCoefficients(); + SH_CONVENTION GetShConvention() const; + void SetShConvention(SH_CONVENTION ShConvention); + protected: ShImage(); ~ShImage() override; template void Construct() const; void PrintSelf(std::ostream &os, itk::Indent indent) const override; mutable mitk::Image::Pointer m_RgbImage; unsigned int m_ShOrder; unsigned int m_NumCoefficients; + SH_CONVENTION m_ShConvention; // use mrtrix style SH convention }; } // namespace mitk #endif /* __mitkShImage__h */ diff --git a/Modules/DiffusionCore/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionCore/mitkDiffusionFunctionCollection.cpp index eaa009d..931a23f 100644 --- a/Modules/DiffusionCore/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionCore/mitkDiffusionFunctionCollection.cpp @@ -1,707 +1,713 @@ /*=================================================================== 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 "mitkDiffusionFunctionCollection.h" #include "mitkNumericTypes.h" #include #include #include #include #include "itkVectorContainer.h" #include "vnl/vnl_vector.h" #include #include #include #include #include #include // Intersect a finite line (with end points p0 and p1) with all of the // cells of a vtkImageData std::vector< std::pair< itk::Index<3>, double > > mitk::imv::IntersectImage(const itk::Vector& spacing, itk::Index<3>& si, itk::Index<3>& ei, itk::ContinuousIndex& sf, itk::ContinuousIndex& ef) { std::vector< std::pair< itk::Index<3>, double > > out; if (si == ei) { double d[3]; for (int i=0; i<3; ++i) d[i] = static_cast(sf[i]-ef[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); out.push_back( std::pair< itk::Index<3>, double >(si, len) ); return out; } double bounds[6]; double entrancePoint[3]; double exitPoint[3]; double startPoint[3]; double endPoint[3]; double t0, t1; for (unsigned int i=0; i<3; ++i) { startPoint[i] = static_cast(sf[i]); endPoint[i] = static_cast(ef[i]); if (si[i]>ei[i]) { auto t = si[i]; si[i] = ei[i]; ei[i] = t; } } for (auto x = si[0]; x<=ei[0]; ++x) for (auto y = si[1]; y<=ei[1]; ++y) for (auto z = si[2]; z<=ei[2]; ++z) { bounds[0] = static_cast(x) - 0.5; bounds[1] = static_cast(x) + 0.5; bounds[2] = static_cast(y) - 0.5; bounds[3] = static_cast(y) + 0.5; bounds[4] = static_cast(z) - 0.5; bounds[5] = static_cast(z) + 0.5; int entryPlane; int exitPlane; int hit = vtkBox::IntersectWithLine(bounds, startPoint, endPoint, t0, t1, entrancePoint, exitPoint, entryPlane, exitPlane); if (hit) { if (entryPlane>=0 && exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i] - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (entryPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (static_cast(ef[i]) - entrancePoint[i])*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } else if (exitPlane>=0) { double d[3]; for (int i=0; i<3; ++i) d[i] = (exitPoint[i]-static_cast(sf[i]))*spacing[i]; double len = std::sqrt( d[0]*d[0] + d[1]*d[1] + d[2]*d[2] ); itk::Index<3> idx; idx[0] = x; idx[1] = y; idx[2] = z; out.push_back( std::pair< itk::Index<3>, double >(idx, len) ); } } } return out; } //------------------------- SH-function ------------------------------------ double mitk::sh::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } void mitk::sh::Cart2Sph(double x, double y, double z, double *spherical) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); if( rad < mitk::eps ) { th = itk::Math::pi/2; phi = itk::Math::pi/2; } else { th = acos(z/rad); phi = atan2(y, x); } spherical[0] = phi; spherical[1] = th; spherical[2] = rad; } vnl_vector_fixed mitk::sh::Sph2Cart(const double& theta, const double& phi, const double& rad) { vnl_vector_fixed dir; dir[0] = rad * sin(theta) * cos(phi); dir[1] = rad * sin(theta) * sin(phi); dir[2] = rad * cos(theta); return dir; } double mitk::sh::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i(::boost::math::spherical_harmonic_r(static_cast(k), -m, theta, phi)); else if (m==0) return static_cast(::boost::math::spherical_harmonic_r(static_cast(k), m, theta, phi)); else return pow(-1.0,m)*sqrt(2.0)*static_cast(::boost::math::spherical_harmonic_i(static_cast(k), m, theta, phi)); } else { double plm = static_cast(::boost::math::legendre_p(k,abs(m),-cos(theta))); double mag = sqrt((2.0*k+1.0)/(4.0*itk::Math::pi)*::boost::math::factorial(k-abs(m))/::boost::math::factorial(k+abs(m)))*plm; if (m>0) return mag*static_cast(cos(m*phi)); else if (m==0) return mag; else return mag*static_cast(sin(-m*phi)); } return 0; } mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromShImage(mitk::Image::Pointer mitkImage) { mitk::ShImage::Pointer mitkShImage = dynamic_cast(mitkImage.GetPointer()); if (mitkShImage.IsNull()) mitkThrow() << "Input image is not a SH image!"; mitk::OdfImage::ItkOdfImageType::Pointer output; switch (mitkShImage->ShOrder()) { case 2: { typedef itk::ShToOdfImageFilter< float, 2 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); + converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 4: { typedef itk::ShToOdfImageFilter< float, 4 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); + converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 6: { typedef itk::ShToOdfImageFilter< float, 6 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); + converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 8: { typedef itk::ShToOdfImageFilter< float, 8 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); + converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 10: { typedef itk::ShToOdfImageFilter< float, 10 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); + converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } case 12: { typedef itk::ShToOdfImageFilter< float, 12 > ShConverterType; typename ShConverterType::InputImageType::Pointer itkvol = ShConverterType::InputImageType::New(); mitk::CastToItkImage(mitkImage, itkvol); typename ShConverterType::Pointer converter = ShConverterType::New(); converter->SetInput(itkvol); + converter->SetToolkit(mitkShImage->GetShConvention()); converter->Update(); output = converter->GetOutput(); break; } default: mitkThrow() << "SH orders higher than 12 are not supported!"; } return output; } mitk::OdfImage::Pointer mitk::convert::GetOdfFromShImage(mitk::Image::Pointer mitkImage) { mitk::OdfImage::Pointer image = mitk::OdfImage::New(); auto img = GetItkOdfFromShImage(mitkImage); image->InitializeByItk( img.GetPointer() ); image->SetVolume( img->GetBufferPointer() ); return image; } mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromTensorImage(mitk::Image::Pointer mitkImage) { typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( GetItkTensorFromTensorImage(mitkImage) ); filter->Update(); return filter->GetOutput(); } mitk::TensorImage::ItkTensorImageType::Pointer mitk::convert::GetItkTensorFromTensorImage(mitk::Image::Pointer mitkImage) { typedef mitk::ImageToItk< mitk::TensorImage::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImage); caster->Update(); return caster->GetOutput(); } mitk::PeakImage::ItkPeakImageType::Pointer mitk::convert::GetItkPeakFromPeakImage(mitk::Image::Pointer mitkImage) { typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImage); caster->SetCopyMemFlag(true); caster->Update(); return caster->GetOutput(); } mitk::OdfImage::Pointer mitk::convert::GetOdfFromTensorImage(mitk::Image::Pointer mitkImage) { mitk::OdfImage::Pointer image = mitk::OdfImage::New(); auto img = GetItkOdfFromTensorImage(mitkImage); image->InitializeByItk( img.GetPointer() ); image->SetVolume( img->GetBufferPointer() ); return image; } mitk::OdfImage::ItkOdfImageType::Pointer mitk::convert::GetItkOdfFromOdfImage(mitk::Image::Pointer mitkImage) { typedef mitk::ImageToItk< mitk::OdfImage::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImage); caster->Update(); return caster->GetOutput(); } vnl_matrix mitk::sh::CalcShBasisForDirections(unsigned int sh_order, vnl_matrix U, bool mrtrix) { vnl_matrix sh_basis = vnl_matrix(U.cols(), (sh_order*sh_order + sh_order + 2)/2 + sh_order ); for(unsigned int i=0; i(sh_order); k+=2) { for(int m=-k; m<=k; m++) { int j = (k*k + k + 2)/2 + m - 1; double phi = U(0,i); double th = U(1,i); sh_basis(i,j) = mitk::sh::Yj(m,k,th,phi, mrtrix); } } } return sh_basis; } unsigned int mitk::sh::ShOrder(int num_coeffs) { int c=3, d=2-2*num_coeffs; int D = c*c-4*d; if (D>0) { int s = (-c+static_cast(sqrt(D)))/2; if (s<0) s = (-c-static_cast(sqrt(D)))/2; return static_cast(s); } else if (D==0) return static_cast(-c/2); return 0; } float mitk::sh::GetValue(const vnl_vector &coefficients, const int &sh_order, const double theta, const double phi, const bool mrtrix) { float val = 0; for(int k=0; k<=sh_order; k+=2) { for(int m=-k; m<=k; m++) { unsigned int j = static_cast((k*k + k + 2)/2 + m - 1); val += coefficients[j] * mitk::sh::Yj(m, k, theta, phi, mrtrix); } } return val; } float mitk::sh::GetValue(const vnl_vector &coefficients, const int &sh_order, const vnl_vector_fixed &dir, const bool mrtrix) { double spherical[3]; mitk::sh::Cart2Sph(dir[0], dir[1], dir[2], spherical); float val = 0; for(int k=0; k<=sh_order; k+=2) { for(int m=-k; m<=k; m++) { int j = (k*k + k + 2)/2 + m - 1; val += coefficients[j] * mitk::sh::Yj(m, k, spherical[1], spherical[0], mrtrix); } } return val; } //------------------------- gradients-function ------------------------------------ mitk::gradients::GradientDirectionContainerType::Pointer mitk::gradients::ReadBvalsBvecs(std::string bvals_file, std::string bvecs_file, double& reference_bval) { mitk::gradients::GradientDirectionContainerType::Pointer directioncontainer = mitk::gradients::GradientDirectionContainerType::New(); std::vector bvec_entries; if (!itksys::SystemTools::FileExists(bvecs_file)) mitkThrow() << "bvecs file not existing: " << bvecs_file; else { std::string line; std::ifstream myfile (bvecs_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; boost::split(strs,line,boost::is_any_of("\t \n")); for (auto token : strs) { if (!token.empty()) { try { bvec_entries.push_back(boost::lexical_cast(token)); } catch(...) { mitkThrow() << "Encountered invalid bvecs file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvecs file could not be opened: " << bvals_file; } } reference_bval = -1; std::vector bval_entries; if (!itksys::SystemTools::FileExists(bvals_file)) mitkThrow() << "bvals file not existing: " << bvals_file; else { std::string line; std::ifstream myfile (bvals_file.c_str()); if (myfile.is_open()) { while (std::getline(myfile, line)) { std::vector strs; boost::split(strs,line,boost::is_any_of("\t \n")); for (auto token : strs) { if (!token.empty()) { try { bval_entries.push_back(boost::lexical_cast(token)); if (bval_entries.back()>reference_bval) reference_bval = bval_entries.back(); } catch(...) { mitkThrow() << "Encountered invalid bvals file entry >" << token << "<"; } } } } myfile.close(); } else { mitkThrow() << "bvals file could not be opened: " << bvals_file; } } for(unsigned int i=0; i0) { double factor = b_val/reference_bval; if(vec.magnitude() > 0) { vec.normalize(); vec[0] = sqrt(factor)*vec[0]; vec[1] = sqrt(factor)*vec[1]; vec[2] = sqrt(factor)*vec[2]; } } directioncontainer->InsertElement(i,vec); } return directioncontainer; } void mitk::gradients::WriteBvalsBvecs(std::string bvals_file, std::string bvecs_file, GradientDirectionContainerType::Pointer gradients, double reference_bval) { std::ofstream myfile; myfile.open (bvals_file.c_str()); for(unsigned int i=0; iSize(); i++) { double twonorm = gradients->ElementAt(i).two_norm(); myfile << std::round(reference_bval*twonorm*twonorm) << " "; } myfile.close(); std::ofstream myfile2; myfile2.open (bvecs_file.c_str()); for(int j=0; j<3; j++) { for(unsigned int i=0; iSize(); i++) { GradientDirectionType direction = gradients->ElementAt(i); direction.normalize(); myfile2 << direction.get(j) << " "; } myfile2 << std::endl; } } std::vector mitk::gradients::GetAllUniqueDirections(const BValueMap & refBValueMap, GradientDirectionContainerType *refGradientsContainer ) { IndiciesVector directioncontainer; auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); auto containerIt = directioncontainer.begin(); bool directionExist = false; while(containerIt != directioncontainer.end()) { if (fabs(dot_product(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { directioncontainer.push_back(wntIndex); } } } return directioncontainer; } bool mitk::gradients::CheckForDifferingShellDirections(const BValueMap & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer) { auto mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ auto mapIterator_2 = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator_2++; //skip bzero Values for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){ if(mapIterator_2 == mapIterator) continue; IndiciesVector currentShell = mapIterator->second; IndiciesVector testShell = mapIterator_2->second; for (unsigned int i = 0; i< currentShell.size(); i++) if (fabs(dot_product(refGradientsContainer->ElementAt(currentShell[i]), refGradientsContainer->ElementAt(testShell[i]))) <= 0.9998) { return true; } } } return false; } vnl_matrix mitk::gradients::ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer) { vnl_matrix Q(3, refShell.size()); Q.fill(0.0); for(unsigned int i = 0; i < refShell.size(); i++) { GradientDirectionType dir = refGradientsContainer->ElementAt(refShell[i]); double x = dir.normalize().get(0); double y = dir.normalize().get(1); double z = dir.normalize().get(2); double cart[3]; mitk::sh::Cart2Sph(x,y,z,cart); Q(0,i) = cart[0]; Q(1,i) = cart[1]; Q(2,i) = cart[2]; } return Q; } vnl_matrix mitk::gradients::ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder) { vnl_matrix SHBasisOutput(QBallReference.cols(), (LOrder+1)*(LOrder+2)*0.5); SHBasisOutput.fill(0.0); for(unsigned int i=0; i 1){ mapIterator++; //skip bzero Values vnl_vector_fixed vec; vec.fill(0.0); directioncontainer->push_back(vec); } for( ; mapIterator != bValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); mitk::gradients::GradientDirectionContainerType::Iterator containerIt = directioncontainer->Begin(); bool directionExist = false; while(containerIt != directioncontainer->End()) { if (fabs(dot_product(containerIt.Value(), origninalGradentcontainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { GradientDirectionType dir(origninalGradentcontainer->ElementAt(wntIndex)); directioncontainer->push_back(dir.normalize()); } } } return directioncontainer; } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkShImageReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkShImageReader.cpp index 83f4dd5..7980963 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkShImageReader.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkShImageReader.cpp @@ -1,131 +1,141 @@ /*=================================================================== 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 "mitkShImageReader.h" #include #include "mitkDiffusionIOMimeTypes.h" #include "itkImageFileReader.h" #include "itkImageRegionIterator.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include "itkNiftiImageIO.h" #include "mitkITKImageImport.h" #include "mitkImageDataItem.h" #include #include #include namespace mitk { ShImageReader::ShImageReader(const ShImageReader& other) : mitk::AbstractFileReader(other) { } ShImageReader::ShImageReader() : mitk::AbstractFileReader( CustomMimeType( mitk::DiffusionIOMimeTypes::SH_MIMETYPE() ), mitk::DiffusionIOMimeTypes::SH_MIMETYPE_DESCRIPTION() ) { + Options defaultOptions; + defaultOptions["Assume MRtrix/MITK style SH convention (else FSL/Dipy)"] = true; + this->SetDefaultOptions(defaultOptions); + m_ServiceReg = this->RegisterService(); } ShImageReader::~ShImageReader() { } template mitk::Image::Pointer ShImageReader::ConvertShImage(ShImage::ShOnDiskType::Pointer img) { typedef itk::ShCoefficientImageImporter< float, shOrder > ImporterType; typename ImporterType::Pointer importer = ImporterType::New(); importer->SetInputImage(img); importer->GenerateData(); mitk::ShImage::Pointer shImage = mitk::ShImage::New(); + + Options options = this->GetOptions(); + bool mrtrix_sh = us::any_cast(options["Assume MRtrix/MITK style SH convention (else FSL/Dipy)"]); + if (!mrtrix_sh) + shImage->SetShConvention(mitk::ShImage::SH_CONVENTION::FSL); + mitk::Image::Pointer resultImage = dynamic_cast(shImage.GetPointer()); mitk::CastToMitkImage(importer->GetCoefficientImage(), resultImage); resultImage->SetVolume(importer->GetCoefficientImage()->GetBufferPointer()); return resultImage; } std::vector > ShImageReader::Read() { mitk::LocaleSwitch localeSwitch("C"); std::vector > result; std::string location = GetInputLocation(); if ( location == "") { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename is empty!"); } else { try { std::string ext = itksys::SystemTools::GetFilenameExtension(location); typedef itk::ImageFileReader< ShImage::ShOnDiskType > FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(location); if (ext==".shi") { itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); reader->SetImageIO(io); } reader->Update(); ShImage::ShOnDiskType::Pointer img = reader->GetOutput(); switch (img->GetLargestPossibleRegion().GetSize()[3]) { case 6: result.push_back( ConvertShImage<2>(img).GetPointer() ); break; case 15: result.push_back( ConvertShImage<4>(img).GetPointer() ); break; case 28: result.push_back( ConvertShImage<6>(img).GetPointer() ); break; case 45: result.push_back( ConvertShImage<8>(img).GetPointer() ); break; case 66: result.push_back( ConvertShImage<10>(img).GetPointer() ); break; case 91: result.push_back( ConvertShImage<12>(img).GetPointer() ); break; default : mitkThrow() << "SH order larger 12 not supported"; } } catch(std::exception& e) { throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what()); } catch(...) { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested vessel tree file!"); } } return result; } } //namespace MITK mitk::ShImageReader* mitk::ShImageReader::Clone() const { return new ShImageReader(*this); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/dipy_reconstructions.py b/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/dipy_reconstructions.py index c9792c2..080abec 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/dipy_reconstructions.py +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.python/resources/dipy_reconstructions.py @@ -1,424 +1,425 @@ import sys - def get_mitk_sphere(): """ Return MITK compliant dipy Sphere object. MITK stores ODFs as 252 values spherically sampled from the continuous ODF. The sampling directions are generate by a 5-fold subdivisions of an icosahedron. """ xyz = np.array([ 0.9756767549555488, 0.9977154378498742, 0.9738192119472443, 0.8915721200771204, 0.7646073555341725, 0.6231965669156312, 0.9817040172417226, 0.9870396762453547, 0.9325589150767597, 0.8173592116492303, 0.6708930871960926, 0.9399233672993689, 0.9144882783890762, 0.8267930935417315, 0.6931818659696647, 0.8407280774774689, 0.782394344826989, 0.6762337155773353, 0.7005607434301688, 0.6228579759074076, 0.5505632701289798, 0.4375940376503738, 0.3153040621970065, 0.1569517536476641, -0.01984099037382634, -0.1857690950088067, -0.3200730131503601, 0.5232435944036425, 0.3889403678268736, 0.2135250052622625, 0.02420694871807206, -0.1448539951504302, 0.5971534158422009, 0.4482053228282282, 0.2597018771197477, 0.06677517278138323, 0.6404616222418184, 0.4782876117785159, 0.2868761951248767, 0.6459894362878276, 0.4789651252338281, 0.3200724178002418, 0.4973180497018747, 0.6793811951363423, 0.8323587928990375, 0.9308933612987835, 0.4036036036586492, 0.5984781165037405, 0.7817280923310203, 0.9140795130247613, 0.4809905907165384, 0.6759621154318279, 0.8390728924802671, 0.5347729120192694, 0.7094340284155564, 0.5560356639783846, 0.2502538949373057, 0.3171352000240629, 0.3793963897789465, 0.4231100429674418, 0.4410301813437042, 0.4357529867703999, 0.5208717223808415, 0.5850086433327374, 0.611055499882272, 0.6009463532173235, 0.6305067000562991, 0.7188806066405239, 0.7654898954879897, 0.7616477696596397, 0.7997756996573342, 0.8700831379830764, 0.8872031228985237, 0.9155019734809123, 0.9568003701205341, -0.4375932291383153, -0.3153035222278598, -0.1569515927579475, 0.0198407706589918, 0.1857686171195431, -0.2644927501381796, -0.1064219080255857, 0.07849995612144045, 0.2583107784678281, -0.04938676750055992, 0.1358448755096817, 0.3243479900672576, 0.1811879481039926, 0.3692668145365748, 0.3890115016151001, -0.6231952788307174, -0.4943551945928708, -0.319458133528771, -0.1156489798772063, 0.08328895892415776, -0.4789641985801549, -0.3127252940830145, -0.1059392282183739, 0.1077444781964869, 0.2912280153186658, -0.2868758523956744, -0.08856892011805101, 0.1287405357080231, 0.3245517154572714, -0.06677541204276306, 0.1413542883070481, 0.3408430926744944, 0.1448534358763926, 0.3374016489097037, -0.2502532651867688, -0.3171345072414974, -0.3793956104585266, -0.4231091882680272, -0.4410293135613324, -0.09929959410007272, -0.1535127609134815, -0.2052877394623771, -0.2436963810571767, 0.08175409117371149, 0.04056025153798869, -0.006048944565669369, 0.2686152102237028, 0.2319923070602857, 0.430309819720559, -0.975676581463901, -0.9977153903038788, -0.9738191090293654, -0.8915716840571059, -0.7646064477130079, -0.9568001079664734, -0.9598482725023617, -0.9044523389503778, -0.7901672201648241, -0.6459882395962464, -0.8872027729137049, -0.8582754834679532, -0.7705800268610806, -0.6404605781008121, -0.7616472974254324, -0.7008201753656432, -0.5971525097007422, -0.6009457148226922, -0.5232427588825813, 0.4943566966479628, 0.3194596781650836, 0.1156503154178581, -0.0832879858164388, 0.5222841738261358, 0.3225497922064885, 0.1018140973507329, 0.5217885230992481, 0.3044789836562512, 0.4873191346491355, -0.4973183240635209, -0.6793811856410323, -0.8323586364840968, -0.9308931819742911, -0.3374020539278631, -0.5261951664998159, -0.7070125356849136, -0.8417962075837926, -0.9155017573317124, -0.3408433114184408, -0.5265312606271311, -0.6896418460594331, -0.7997755164970677, -0.3245517106425898, -0.4925847482169691, -0.6305065080228541, -0.2912277152063287, -0.4357526334612896, 0.7901679726328494, 0.9044526665335126, 0.9598484396937114, 0.7705806468939737, 0.858275831469383, 0.7008207681995118, -0.4036039458806759, -0.2583110138480089, -0.0784999126587471, 0.1064223584250461, 0.264493571710179, -0.4809907334514471, -0.3243480295764106, -0.1358446002697818, 0.04938746901646566, -0.5347730026038946, -0.3692667658371347, -0.1811875286592425, -0.5560358190148772, -0.3890114324926668, -0.5505634949474449, 0.8417963565884857, 0.7070125813068046, 0.5261950179989611, 0.6896418985458221, 0.5265311900255359, 0.4925848265160583, 0.2436972866599269, 0.2052886581368649, 0.153513629451971, 0.09930039009433847, 0.006049691633511915, -0.04055950638179381, -0.08175337578691833, -0.2319919155781195, -0.2686148310916902, -0.430309819678344, -0.02420720081803753, -0.2135248270679241, -0.3889397838050994, -0.2597016312374675, -0.4482046405142344, -0.4782867918076852, -0.1018130528605821, -0.322548598821141, -0.5222830294256716, -0.6708921376896406, -0.304478224282928, -0.5217878437313506, -0.6931813485878851, -0.4873188675145023, -0.6762335873429084, -0.6228580878699612, -0.6110548409057, -0.5850080622199078, -0.5208712693637837, -0.7654894328832393, -0.7188802647693375, -0.8700828159137221, -0.8173587433845655, -0.9325588839421305, -0.9870397834787261, -0.9817039872478999, -0.8267930492778305, -0.9144884914916022, -0.9399235077793813, -0.7823945479956939, -0.8407283372889187, -0.7005610213599369, -0.1077438933887955, 0.1059400956623477, 0.3127262866621893, -0.1287403742204129, 0.08856921814263634, -0.1413545191115968, -0.9140794058749131, -0.7817279594934516, -0.5984781448346268, -0.8390728949381593, -0.6759620794963979, -0.709434131000089, -0.1778161375899129, -0.06053925384414331, 0.07929679392711581, 0.222673458561735, 0.3458247516791153, 0.4366423972091846, 0.01030826616734189, 0.1591522280204451, 0.3173816763430465, 0.4549463955350546, 0.5521270265729551, 0.2292788658415479, 0.3973400932411465, 0.5502139834879405, 0.6594089221868847, 0.4476465561008348, 0.6096570464011057, 0.7343998566036512, 0.629214796874201, 0.7646693979379596, 0.7580253719184178, -0.5980610514568761, -0.5101530988159087, -0.382225667160838, -0.2244621267538426, -0.06301328229424107, 0.07805400320688782, -0.4311039309963852, -0.3079662136138592, -0.1501157132113724, 0.01750888497279251, 0.1650825345160538, -0.2148810450151756, -0.06090095222676627, 0.1073128739652992, 0.2584097661066967, 0.02655484252908358, 0.1901297170957776, 0.3420822257932489, 0.2531835106264871, 0.4022303494272352, -0.07805410188638827, -0.1080255529483224, -0.1376217050758367, -0.1609000070073124, -0.1740018618448228, 0.09827676798573926, 0.083291898217249, 0.06127443921955168, 0.03526739273256396, 0.2991139104294396, 0.2941068360088736, 0.2692865316145088, 0.4942032775296958, 0.4857723178878524, 0.6512069539966677, -0.9161616153729886, -0.9396953110011561, -0.9204280785344878, -0.8462030522374957, -0.7293237120999879, -0.8470541513588044, -0.8482966176587544, -0.7977006542517769, -0.6951661565374421, -0.566558592627622, -0.7243096319272092, -0.6931460376496088, -0.6140043047773551, -0.5016343691560573, -0.5520254073275178, -0.4928644880867128, -0.403575153350467, -0.3587591578566765, -0.2886351685087218, 0.5980613131647216, 0.5101532951859686, 0.382225843595672, 0.2244622808787926, 0.06301334452030186, 0.6944632949786616, 0.5955168212825119, 0.4473425940100297, 0.2700417838303327, 0.7724043956082883, 0.6553545192922715, 0.4871408620353512, 0.8097301284690857, 0.6725220182496192, 0.8002534097038426, -0.4366431953633789, -0.5869882376922511, -0.7332080507197046, -0.8450980113065225, -0.9041113586460733, -0.4022310083998925, -0.554596445154436, -0.6925605687496104, -0.7854318984598006, -0.8250621271173465, -0.3420827953668352, -0.4840440064641756, -0.6033456975789954, -0.6777531805937266, -0.2584102557043402, -0.3819753792546441, -0.4821906665520286, -0.1650828712784331, -0.270790845781693, 0.9161619708326184, 0.9396956424389374, 0.9204283182965946, 0.8462032095340455, 0.7293238793541417, 0.9749588444840027, 0.9879501207294071, 0.942053498973333, 0.8348196077814718, 0.9950795014807369, 0.9818515654328379, 0.9027098746674149, 0.9581801446138297, 0.9118246030313639, 0.8703772282258925, 0.1778171059718252, 0.06053992567271226, -0.07929659020903117, -0.2226737578340799, -0.345825401239635, 0.2886360377097776, 0.1672516508448342, 0.02000533874392893, -0.1285435155191929, -0.2531843553864728, 0.403575906447316, 0.2774342678683828, 0.1245598363284875, -0.02655554762561945, 0.5016349858535857, 0.3695530582277636, 0.2148806720954671, 0.5665590425344393, 0.431103930292903, 0.5869876102086139, 0.7332077514676827, 0.845098078457225, 0.9041116580482536, 0.7182616282077119, 0.8617334421407644, 0.9490975365686583, 0.8223898048944452, 0.9416915744235097, 0.8729720010540123, 0.1080256414522809, 0.1376220280275969, 0.1609005865750696, 0.1740026689030255, 0.2707904196202965, 0.3196768235430837, 0.3552546724685221, 0.3677018240803483, 0.3587598208776521, 0.4821901792282771, 0.5389508449256169, 0.5637713635689835, 0.5520258363563475, 0.6777529577987501, 0.7231337276202411, 0.724309982145211, 0.8250622687013296, 0.8470545173149734, 0.1285429999155006, -0.02000532948058562, -0.1672511147059996, -0.1245600244829796, -0.2774338902981233, -0.3695528631494325, -0.09827641615811868, -0.2700412859530667, -0.4473420975374328, -0.5955164071695848, -0.6944629164413806, -0.2991130971968019, -0.4871400501186961, -0.6553538941234454, -0.7724039524031648, -0.4942022299541438, -0.6725212074710563, -0.8097296395344389, -0.6512059956089504, -0.8002528392148971, -0.7580246814516085, -0.3677014077761052, -0.3552545716101517, -0.3196770257819652, -0.5637712030900536, -0.5389510214534028, -0.7231336172569296, -0.8348194119106425, -0.9420533966954356, -0.9879499956150448, -0.9749586635216289, -0.9027097279159257, -0.9818515951566739, -0.9950795477220543, -0.9118244750171576, -0.9581802235578871, -0.8703770126934449, -0.0175091339170676, 0.1501155140512474, 0.3079660822386824, -0.1073133727582037, 0.06090046334304851, -0.1901304002696938, -0.9490974969653682, -0.8617336589899791, -0.7182621005240754, -0.5521276321758419, -0.941691783045487, -0.8223901593137167, -0.6594093292610237, -0.872972144171723, -0.7343999908188845, -0.7646691446910742, 0.6951665021597787, 0.7977009700656229, 0.8482969664746548, 0.6140047811934269, 0.6931464276818936, 0.4928650597255946, -0.4549467775084718, -0.3173815862988101, -0.1591515620353438, -0.01030716362341688, -0.5502140363721867, -0.3973395475484636, -0.2292777334167206, -0.609656670182737, -0.4476455277450017, -0.6292139442700462, 0.7854319364049284, 0.6925603649758249, 0.5545959620739339, 0.6033453603619342, 0.4840435291285519, 0.3819748711371402, -0.03526641653115874, -0.06127364342066123, -0.0832913202753871, -0.2692854573778917, -0.2941058574917593, -0.4857712605383084, -0.1282040992012934, 0.02998172476739921, 0.2130449739264662, 0.394354771181159, 0.5438573645627299, 0.6488215902264565, -0.1901340637026595, -0.02057293935230464, 0.1720544722828635, 0.3534794142829396, 0.4950335464190314, -0.252933321812349, -0.07636778766011496, 0.1169519253626288, 0.2909961752861106, -0.304612640171253, -0.1271903099934383, 0.0580070042064605, -0.3366056805211806, -0.1653138037361849, -0.3496821715684524, 0.6714420577705413, 0.8002044514563711, 0.9106424580428781, 0.9742808059046055, 0.9805708386415104, 0.9441720387917811, 0.7350956003099328, 0.8682639008659977, 0.9653353535299492, 0.9995536316680411, 0.9755844796257857, 0.7728091190204586, 0.8918537226509272, 0.9597077065970592, 0.9637247890765801, 0.767530944505584, 0.857374860312736, 0.8948082473172733, 0.7201359303293944, 0.7802583897718675, -0.9441722324566724, -0.8608166107545396, -0.7207644955095487, -0.5303678229575245, -0.3211867088850157, -0.9096404828216634, -0.7967975927156801, -0.620601831095295, -0.4039985827676406, -0.8241231220089414, -0.6757043639889994, -0.4726959329165278, -0.6854057579633669, -0.5106159168102177, -0.5164821811548767, -0.3130828685601688, -0.128054626578418, 0.09418349997750548, 0.3239109228229815, 0.523048087763098, -0.3043330399192493, -0.09531787499064681, 0.146419102006115, 0.3786227553496849, 0.5638039035645359, -0.2789925774668332, -0.05252850546893189, 0.1924160430438771, 0.4101897544477446, -0.2358533016570041, -0.006318969895916147, 0.2236016867495729, -0.1820659309330632, 0.03496843200875603, -0.6714423515894649, -0.8002045390283683, -0.9106424117173109, -0.9742807748705438, -0.9805709251787237, -0.6691519386893557, -0.79626257796142, -0.8909109722488041, -0.9275521423149625, -0.6332080201962388, -0.7430051304270751, -0.8108589038018792, -0.5581290590099237, -0.6413705283620924, -0.4563600901355626, -0.648822290296783, -0.6411378559950974, -0.600293640880965, -0.5219527418637842, -0.4191009430775375, -0.7802586188950914, -0.7711197529973893, -0.713538182957094, -0.6094980396194888, -0.4842093859996422, -0.8948081394501657, -0.8705497953564927, -0.7870195954857328, -0.6597854273844109, -0.9637246412193355, -0.9132982945459158, -0.8070428410352181, -0.975584505681221, -0.9015722073987464, 0.3130823317650696, 0.1280539101236855, -0.09418429615654819, -0.3239116283455282, -0.523048586255087, 0.1989845274738189, -0.01970764286946627, -0.2653151882168217, -0.4936479477555078, 0.05597369301027853, -0.1852629737758222, -0.4302072668465533, -0.09867461327158224, -0.3387557568094476, -0.2393260112020272, 0.1282040764044601, -0.0299819504088954, -0.2130455201840074, -0.3943555879655132, -0.5438582278251758, -0.03496843048443223, -0.2252069693946209, -0.4261053308619027, -0.5992598174372795, -0.720136706807033, -0.2236017161594696, -0.4317330442416767, -0.6250530132529536, -0.7675317913865697, -0.4101898771205939, -0.6101488498350025, -0.7728099228904255, -0.5638041319099237, -0.7350961954485662, 0.6411372723067146, 0.6002931843810706, 0.5219523372221417, 0.4191004905746099, 0.4596949806286311, 0.3916338931244087, 0.2980734064957148, 0.226741584116328, 0.1432114770939381, -0.02097489882147943, 0.8608164411414747, 0.7207644427956729, 0.5303678926086439, 0.3211867913977836, 0.9015721838250796, 0.7879881821713033, 0.6114960278478284, 0.3951892122332402, 0.1820657113417612, 0.8070430398170311, 0.6574928275930984, 0.4544842943197335, 0.2358529185889448, 0.6597856586149884, 0.4841878538357612, 0.2789921022280572, 0.4842093252521232, 0.3043325272261384, 0.5992589358516202, 0.4261046359672609, 0.2252066549797059, 0.6250522113657903, 0.4317325950511361, 0.6101482870567641, 0.9096403689902206, 0.9275522217134882, 0.8909112253661301, 0.796262827475376, 0.6691520068054228, 0.8241233338640371, 0.810859375773786, 0.7430057321681839, 0.6332085061147845, 0.6854064426268304, 0.6413714065577412, 0.5581299045184589, 0.5164832226272315, 0.4563611494403301, 0.3496833143594963, -0.3951892821849063, -0.6114960336943951, -0.787988199289983, -0.4544844137443082, -0.657492739431111, -0.484187939006181, 0.4936478319326018, 0.2653148405479006, 0.01970714938077021, -0.1989850169013517, 0.4302075642722875, 0.1852629793843341, -0.0559739158243807, 0.3387563694841473, 0.09867487876932232, 0.2393267951217032, -0.999553621201999, -0.9653354239158236, -0.8682642090770526, -0.9597077173597477, -0.8918540989344099, -0.8573751662344773, -0.2980738893651726, -0.3916343988495664, -0.4596955428592778, -0.4950341577852201, -0.1432117197792371, -0.2267418620329016, -0.2909964852939082, 0.02097514873862574, -0.05800679989935065, 0.1653145532988453, -0.3786231842883476, -0.1464197032303796, 0.09531724619007391, -0.1924163631703616, 0.05252803743712917, 0.006318730357784829, -0.3534800054422614, -0.1720548071373146, 0.02057294660420643, 0.190134278339324, -0.1169519894866824, 0.07636807502743861, 0.2529338262925594, 0.1271908635410245, 0.3046134343217798, 0.3366066958443542, 0.6094980941008995, 0.7135382519498201, 0.7711196978950583, 0.7870198804193677, 0.8705500304441893, 0.9132984713369965, 0.403998910419839, 0.62060207699311, 0.7967976318501995, 0.4726965405256068, 0.6757048258462731, 0.5106167801856609]) n = int(xyz.shape[0] / 3) x = xyz[:n] y = xyz[n:2 * n] z = xyz[2 * n:] for i in range(n): v = np.array([x[i], y[i], z[i]]) norm = np.linalg.norm(v) if norm > 0: v /= norm x[i] = v[0] y[i] = v[1] z[i] = v[2] s = sphere.Sphere(x=x, y=y, z=z) return s error_string = None del error_string + try: import dipy.direction.peaks as dpp from dipy.reconst.shore import ShoreModel from dipy.reconst.shm import CsaOdfModel, OpdtModel, SphHarmModel import dipy.reconst.sfm as sfm from dipy.reconst.csdeconv import auto_response, ConstrainedSphericalDeconvModel from dipy.core import sphere import numpy as np from dipy.core.gradients import gradient_table import SimpleITK as sitk print('DIPY Reconstructions') data = sitk.GetArrayFromImage(in_image) + print(data.shape) bvals = np.array(bvals) bvecs = np.array(bvecs) # create dipy Sphere sphere = get_mitk_sphere() odf = None model = None gtab = gradient_table(bvals, bvecs) if mask is not None: mask = sitk.GetArrayFromImage(mask) print(mask.shape) # fit selected model sh_coeffs = None odf = None if model_type == '3D-SHORE': print('Fitting 3D-SHORE') print("radial_order: ", radial_order) print("zeta: ", zeta) print("lambdaN: ", lambdaN) print("lambdaL: ", lambdaL) model = ShoreModel(gtab, radial_order=radial_order, zeta=zeta, lambdaN=lambdaN, lambdaL=lambdaL) asmfit = model.fit(data) odf = asmfit.odf(sphere) elif model_type == 'CSA-QBALL': print('Fitting CSA-QBALL') print("sh_order: ", sh_order) print("smooth: ", smooth) model = CsaOdfModel(gtab=gtab, sh_order=sh_order, smooth=smooth) sh_coeffs = model.fit(data, mask=mask).shm_coeff elif model_type == 'SFM': print('Fitting SFM') print("fa_thr: ", fa_thr) response, ratio = auto_response(gtab, data, roi_radius=10, fa_thr=fa_thr) model = sfm.SparseFascicleModel(gtab, sphere=sphere, l1_ratio=0.5, alpha=0.001, response=response[0]) odf = model.fit(data, mask=mask).odf(sphere) elif model_type == 'CSD': print('Fitting CSD') print("sh_order: ", sh_order) print("fa_thr: ", fa_thr) response, ratio = auto_response(gtab, data, roi_radius=10, fa_thr=fa_thr) model = ConstrainedSphericalDeconvModel(gtab, response, sh_order=sh_order) sh_coeffs = model.fit(data).shm_coeff elif model_type == 'Opdt': print('Orientation Probability Density Transform') print("sh_order: ", sh_order) print("smooth: ", smooth) model = OpdtModel(gtab=gtab, sh_order=sh_order, smooth=smooth) sh_coeffs = model.fit(data, mask=mask).shm_coeff else: raise ValueError('Model type not supported. Available models: 3D-SHORE, CSA-QBALL, SFM, CSD, Opdt') if odf is not None: odf = np.nan_to_num(odf) print('Preparing ODF image') odf_image = sitk.Image([data.shape[2], data.shape[1], data.shape[0]], sitk.sitkVectorFloat32, len(sphere.vertices)) for x in range(data.shape[2]): for y in range(data.shape[1]): for z in range(data.shape[0]): if mask is not None and mask[z, y, x] == 0: continue odf_image.SetPixel(x, y, z, odf[z, y, x, :]) odf_image.SetOrigin(in_image.GetOrigin()) odf_image.SetSpacing(in_image.GetSpacing()) odf_image.SetDirection(in_image.GetDirection()) elif sh_coeffs is not None: sh_coeffs = np.nan_to_num(sh_coeffs) print('Preparing SH image') sh_image = sitk.Image([sh_coeffs.shape[2], sh_coeffs.shape[1], sh_coeffs.shape[0]], sitk.sitkVectorFloat32, sh_coeffs.shape[3]) for x in range(sh_coeffs.shape[2]): for y in range(sh_coeffs.shape[1]): for z in range(sh_coeffs.shape[0]): if mask is not None and mask[z, y, x] == 0: continue sh_image.SetPixel(x, y, z, sh_coeffs[z, y, x, :]) sh_image.SetOrigin(in_image.GetOrigin()) sh_image.SetSpacing(in_image.GetSpacing()) sh_image.SetDirection(in_image.GetDirection()) if num_peaks > 0: print('Calculating peaks') sys.stdout.flush() # calculate peak image data = np.nan_to_num(data) sf_peaks = dpp.peaks_from_model(model, data, sphere, relative_peak_threshold=relative_peak_threshold, min_separation_angle=min_separation_angle, return_sh=False, npeaks=num_peaks, parallel=True, mask=mask) # reshape to be MITK/MRtrix compliant s = sf_peaks.peak_dirs.shape peaks = sf_peaks.peak_dirs.reshape((s[0], s[1], s[2], num_peaks * 3), order='C') peaks = np.nan_to_num(peaks) peak_image = sitk.Image([data.shape[2], data.shape[1], data.shape[0]], sitk.sitkVectorFloat32, num_peaks * 3) # scale peaks max_peak = 1.0 if normalize_peaks: max_peak = np.max(sf_peaks.peak_values) if max_peak <= 0: max_peak = 1.0 for x in range(s[0]): for y in range(s[1]): for z in range(s[2]): for i in range(num_peaks): peaks[x, y, z, i * 3] *= sf_peaks.peak_values[x, y, z, i] / max_peak peaks[x, y, z, i * 3 + 1] *= sf_peaks.peak_values[x, y, z, i] / max_peak peaks[x, y, z, i * 3 + 2] *= sf_peaks.peak_values[x, y, z, i] / max_peak peak_image.SetPixel(z, y, x, peaks[x, y, z, :]) peak_image.SetOrigin(in_image.GetOrigin()) peak_image.SetSpacing(in_image.GetSpacing()) peak_image.SetDirection(in_image.GetDirection()) except Exception as e: error_string = str(e) print(error_string) sys.stdout.flush() diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.python/src/internal/QmitkDipyReconstructionsView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.python/src/internal/QmitkDipyReconstructionsView.cpp index 9b639bf..c27d339 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.python/src/internal/QmitkDipyReconstructionsView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.python/src/internal/QmitkDipyReconstructionsView.cpp @@ -1,409 +1,410 @@ /*=================================================================== 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 "QmitkDipyReconstructionsView.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkDipyReconstructionsView::VIEW_ID = "org.mitk.views.dipyreconstruction"; QmitkDipyReconstructionsView::QmitkDipyReconstructionsView() : QmitkAbstractView() , m_Controls( 0 ) { } // Destructor QmitkDipyReconstructionsView::~QmitkDipyReconstructionsView() { } void QmitkDipyReconstructionsView::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::QmitkDipyReconstructionsViewControls; m_Controls->setupUi( parent ); connect( m_Controls->m_ImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGUI()) ); connect( m_Controls->m_StartButton, SIGNAL(clicked()), this, SLOT(StartFit()) ); connect( m_Controls->m_ModelBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGUI()) ); this->m_Parent = parent; m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New(); m_Controls->m_ImageBox->SetPredicate( isDwi ); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateAnd::Pointer isBinary3dImage = mitk::NodePredicateAnd::New( mitk::NodePredicateAnd::New( isImagePredicate, isBinaryPredicate ), mitk::NodePredicateDimension::New(3)); m_Controls->m_MaskBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskBox->SetPredicate( isBinary3dImage ); m_Controls->m_MaskBox->SetZeroEntryText("--"); UpdateGUI(); } } void QmitkDipyReconstructionsView::OnSelectionChanged(berry::IWorkbenchPart::Pointer, const QList& ) { } void QmitkDipyReconstructionsView::UpdateGUI() { if (m_Controls->m_ImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_StartButton->setEnabled(true); else m_Controls->m_StartButton->setEnabled(false); m_Controls->m_ShoreBox->setVisible(false); m_Controls->m_SfmBox->setVisible(false); m_Controls->m_CsdBox->setVisible(false); m_Controls->m_CsaBox->setVisible(false); m_Controls->m_OpdtBox->setVisible(false); switch(m_Controls->m_ModelBox->currentIndex()) { case 0: { m_Controls->m_ShoreBox->setVisible(true); break; } case 1: { m_Controls->m_SfmBox->setVisible(true); break; } case 2: { m_Controls->m_CsdBox->setVisible(true); break; } case 3: { m_Controls->m_CsaBox->setVisible(true); break; } case 4: { m_Controls->m_OpdtBox->setVisible(true); break; } } } void QmitkDipyReconstructionsView::SetFocus() { UpdateGUI(); m_Controls->m_StartButton->setFocus(); } void QmitkDipyReconstructionsView::StartFit() { mitk::DataNode::Pointer node = m_Controls->m_ImageBox->GetSelectedNode(); mitk::Image::Pointer input_image = dynamic_cast(node->GetData()); // get python script as string QString data; - QString fileName(":/QmitkDiffusionImaging/dipy_reconstructions.py"); + QString fileName(":/QmitkDiffusionPython/dipy_reconstructions.py"); QFile file(fileName); if(!file.open(QIODevice::ReadOnly)) { - qDebug()<<"filenot opened"<Size(); ++i) { bvals += boost::lexical_cast(bvaluevector.at(i)); if (bvaluevector.at(i)==0) b0_count++; auto g = gcont->GetElement(i); if (g.two_norm()>0.000001) g /= g.two_norm(); bvecs += "[" + boost::lexical_cast(g[0]) + "," + boost::lexical_cast(g[1]) + "," + boost::lexical_cast(g[2]) + "]"; if (iSize()-1) { bvals += ", "; bvecs += ", "; } } bvals += "]"; bvecs += "]"; if (b0_count==0) { QMessageBox::warning(nullptr, "Error", "No b=0 volume found. Do your b-values need rounding? Use the Preprocessing View for rounding b-values,", QMessageBox::Ok); return; } - us::ModuleContext* context = us::GetModuleContext(); us::ServiceReference m_PythonServiceRef = context->GetServiceReference(); mitk::IPythonService* m_PythonService = dynamic_cast ( context->GetService(m_PythonServiceRef) ); mitk::IPythonService::ForceLoadModule(); m_PythonService->CopyToPythonAsSimpleItkImage( input_image, "in_image"); m_PythonService->Execute("mask=None"); if (m_Controls->m_MaskBox->GetSelectedNode().IsNotNull()) { auto mitk_mask = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); if (mitk_mask->GetLargestPossibleRegion().GetSize()==input_image->GetLargestPossibleRegion().GetSize()) m_PythonService->CopyToPythonAsSimpleItkImage( mitk_mask, "mask"); else MITK_INFO << "Mask image not used. Does not match data size: " << mitk_mask->GetLargestPossibleRegion().GetSize() << " vs. " << input_image->GetLargestPossibleRegion().GetSize(); } m_PythonService->Execute("normalize_peaks=False"); if (m_Controls->m_NormalizePeaks->isChecked()) m_PythonService->Execute("normalize_peaks=True"); std::string model = "3D-SHORE"; int sh_order = 0; switch(m_Controls->m_ModelBox->currentIndex()) { case 0: { model = "3D-SHORE"; m_PythonService->Execute("radial_order=" + boost::lexical_cast(m_Controls->m_RadialOrder->value())); m_PythonService->Execute("zeta=" + boost::lexical_cast(m_Controls->m_Zeta->value())); m_PythonService->Execute("lambdaN=" + m_Controls->m_LambdaN->text().toStdString()); m_PythonService->Execute("lambdaL=" + m_Controls->m_LambdaL->text().toStdString()); break; } case 1: { model = "SFM"; m_PythonService->Execute("fa_thr=" + boost::lexical_cast(m_Controls->m_FaThresholdSfm->value())); break; } case 2: { model = "CSD"; sh_order = m_Controls->m_ShOrderCsd->value(); m_PythonService->Execute("sh_order=" + boost::lexical_cast(sh_order)); m_PythonService->Execute("fa_thr=" + boost::lexical_cast(m_Controls->m_FaThresholdCsd->value())); break; } case 3: { model = "CSA-QBALL"; sh_order = m_Controls->m_ShOrderCsa->value(); m_PythonService->Execute("sh_order=" + boost::lexical_cast(sh_order)); m_PythonService->Execute("smooth=" + boost::lexical_cast(m_Controls->m_LambdaCsa->value())); break; } case 4: { model = "Opdt"; sh_order = m_Controls->m_ShOrderOpdt->value(); m_PythonService->Execute("sh_order=" + boost::lexical_cast(sh_order)); m_PythonService->Execute("smooth=" + boost::lexical_cast(m_Controls->m_LambdaOpdt->value())); break; } } m_PythonService->Execute("model_type='"+model+"'"); m_PythonService->Execute("num_peaks=0"); if (m_Controls->m_DoCalculatePeaks->isChecked()) { m_PythonService->Execute("num_peaks=" + boost::lexical_cast(m_Controls->m_NumPeaks->value())); m_PythonService->Execute("min_separation_angle=" + boost::lexical_cast(m_Controls->m_SepAngle->value())); m_PythonService->Execute("relative_peak_threshold=" + boost::lexical_cast(m_Controls->m_RelativeThreshold->value())); } m_PythonService->Execute("data=False"); m_PythonService->Execute("bvals=" + bvals); m_PythonService->Execute("bvecs=" + bvecs); + m_PythonService->Execute(data.toStdString(), mitk::IPythonService::MULTI_LINE_COMMAND); // clean up after running script (better way than deleting individual variables?) if(m_PythonService->DoesVariableExist("in_image")) m_PythonService->Execute("del in_image"); // check for errors if(!m_PythonService->GetVariable("error_string").empty()) { QMessageBox::warning(nullptr, "Error", QString(m_PythonService->GetVariable("error_string").c_str()), QMessageBox::Ok); return; } if (m_PythonService->DoesVariableExist("odf_image")) { mitk::OdfImage::ItkOdfImageType::Pointer itkImg = mitk::OdfImage::ItkOdfImageType::New(); mitk::Image::Pointer out_image = m_PythonService->CopySimpleItkImageFromPython("odf_image"); mitk::CastToItkImage(out_image, itkImg); mitk::OdfImage::Pointer image = mitk::OdfImage::New(); image->InitializeByItk( itkImg.GetPointer() ); image->SetVolume( itkImg->GetBufferPointer() ); mitk::DataNode::Pointer odfs = mitk::DataNode::New(); odfs->SetData( image ); QString name(node->GetName().c_str()); odfs->SetName(name.toStdString() + "_" + model); GetDataStorage()->Add(odfs, node); m_PythonService->Execute("del odf_image"); } else if (m_PythonService->DoesVariableExist("sh_image")) { mitk::Image::Pointer out_image = m_PythonService->CopySimpleItkImageFromPython("sh_image"); itk::VectorImage::Pointer vectorImage = itk::VectorImage::New(); mitk::CastToItkImage(out_image, vectorImage); itk::VectorImageToFourDImageFilter< float >::Pointer converter = itk::VectorImageToFourDImageFilter< float >::New(); converter->SetInputImage(vectorImage); converter->GenerateData(); mitk::ShImage::ShOnDiskType::Pointer itkImg = converter->GetOutputImage(); mitk::ShImage::Pointer shImage = mitk::ShImage::New(); + shImage->SetShConvention(mitk::ShImage::SH_CONVENTION::FSL); mitk::Image::Pointer mitkImage = dynamic_cast(shImage.GetPointer()); switch(sh_order) { case 2: { typedef itk::ShCoefficientImageImporter< float, 2 > ImporterType; typename ImporterType::Pointer importer = ImporterType::New(); importer->SetInputImage(itkImg); importer->GenerateData(); mitk::CastToMitkImage(importer->GetCoefficientImage(), mitkImage); mitkImage->SetVolume(importer->GetCoefficientImage()->GetBufferPointer()); break; } case 4: { typedef itk::ShCoefficientImageImporter< float, 4 > ImporterType; typename ImporterType::Pointer importer = ImporterType::New(); importer->SetInputImage(itkImg); importer->GenerateData(); mitk::CastToMitkImage(importer->GetCoefficientImage(), mitkImage); mitkImage->SetVolume(importer->GetCoefficientImage()->GetBufferPointer()); break; } case 6: { typedef itk::ShCoefficientImageImporter< float, 6 > ImporterType; typename ImporterType::Pointer importer = ImporterType::New(); importer->SetInputImage(itkImg); importer->GenerateData(); mitk::CastToMitkImage(importer->GetCoefficientImage(), mitkImage); mitkImage->SetVolume(importer->GetCoefficientImage()->GetBufferPointer()); break; } case 8: { typedef itk::ShCoefficientImageImporter< float, 8 > ImporterType; typename ImporterType::Pointer importer = ImporterType::New(); importer->SetInputImage(itkImg); importer->GenerateData(); mitk::CastToMitkImage(importer->GetCoefficientImage(), mitkImage); mitkImage->SetVolume(importer->GetCoefficientImage()->GetBufferPointer()); break; } case 10: { typedef itk::ShCoefficientImageImporter< float, 10 > ImporterType; typename ImporterType::Pointer importer = ImporterType::New(); importer->SetInputImage(itkImg); importer->GenerateData(); mitk::CastToMitkImage(importer->GetCoefficientImage(), mitkImage); mitkImage->SetVolume(importer->GetCoefficientImage()->GetBufferPointer()); break; } case 12: { typedef itk::ShCoefficientImageImporter< float, 12 > ImporterType; typename ImporterType::Pointer importer = ImporterType::New(); importer->SetInputImage(itkImg); importer->GenerateData(); mitk::CastToMitkImage(importer->GetCoefficientImage(), mitkImage); mitkImage->SetVolume(importer->GetCoefficientImage()->GetBufferPointer()); break; } default: mitkThrow() << "SH order not supported"; } mitk::DataNode::Pointer shNode = mitk::DataNode::New(); shNode->SetData( mitkImage ); QString name(node->GetName().c_str()); shNode->SetName(name.toStdString() + "_" + model); GetDataStorage()->Add(shNode, node); m_PythonService->Execute("del sh_image"); } if (m_Controls->m_DoCalculatePeaks->isChecked() && m_PythonService->DoesVariableExist("peak_image")) { mitk::Image::Pointer out_image = m_PythonService->CopySimpleItkImageFromPython("peak_image"); itk::VectorImage::Pointer vectorImage = itk::VectorImage::New(); mitk::CastToItkImage(out_image, vectorImage); itk::VectorImageToFourDImageFilter< float >::Pointer converter = itk::VectorImageToFourDImageFilter< float >::New(); converter->SetInputImage(vectorImage); converter->GenerateData(); mitk::PeakImage::ItkPeakImageType::Pointer itk_peaks = converter->GetOutputImage(); mitk::Image::Pointer mitk_peaks = dynamic_cast(mitk::PeakImage::New().GetPointer()); mitk::CastToMitkImage(itk_peaks, mitk_peaks); mitk_peaks->SetVolume(itk_peaks->GetBufferPointer()); mitk::DataNode::Pointer seg = mitk::DataNode::New(); seg->SetData( mitk_peaks ); seg->SetName("Peaks"); GetDataStorage()->Add(seg, node); m_PythonService->Execute("del peak_image"); } + mitk::IPythonService::ForceLoadModule(); // just included here because it flushes the python stdout }