diff --git a/Modules/Connectomics/Rendering/mitkConnectomicsEnumerationSubclassesSerializer.cpp b/Modules/Connectomics/Rendering/mitkConnectomicsEnumerationSubclassesSerializer.cpp index 84d4bd2..b0c6b82 100644 --- a/Modules/Connectomics/Rendering/mitkConnectomicsEnumerationSubclassesSerializer.cpp +++ b/Modules/Connectomics/Rendering/mitkConnectomicsEnumerationSubclassesSerializer.cpp @@ -1,79 +1,80 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkConnectomicsEnumerationSubclassesSerializer_h_included #define mitkConnectomicsEnumerationSubclassesSerializer_h_included #include "mitkEnumerationPropertySerializer.h" #include "mitkConnectomicsRenderingEdgeColorParameterProperty.h" #include "mitkConnectomicsRenderingEdgeFilteringProperty.h" #include "mitkConnectomicsRenderingEdgeRadiusParameterProperty.h" #include "mitkConnectomicsRenderingEdgeThresholdParameterProperty.h" #include "mitkConnectomicsRenderingNodeColorParameterProperty.h" #include "mitkConnectomicsRenderingNodeFilteringProperty.h" #include "mitkConnectomicsRenderingNodeRadiusParameterProperty.h" #include "mitkConnectomicsRenderingNodeThresholdParameterProperty.h" #include "mitkConnectomicsRenderingSchemeProperty.h" +#include #include #define MITK_REGISTER_ENUM_SUB_SERIALIZER(classname) \ \ namespace mitk \ { \ \ class MITKCONNECTOMICS_EXPORT classname ## Serializer : public EnumerationPropertySerializer \ { \ public: \ \ mitkClassMacro( classname ## Serializer, EnumerationPropertySerializer ) \ itkFactorylessNewMacro(Self) \ itkCloneMacro(Self) \ \ - virtual BaseProperty::Pointer Deserialize(TiXmlElement* element) \ + virtual BaseProperty::Pointer Deserialize(const tinyxml2::XMLElement* element) \ { \ if (!element) return nullptr; \ const char* sa( element->Attribute("value") ); \ \ std::string s(sa?sa:""); \ classname::Pointer property = classname::New(); \ property->SetValue( s ); \ \ return property.GetPointer(); \ } \ \ protected: \ \ classname ## Serializer () {} \ virtual ~classname ## Serializer () {} \ }; \ \ } \ \ MITK_REGISTER_SERIALIZER( classname ## Serializer ); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingEdgeColorParameterProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingEdgeFilteringProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingEdgeRadiusParameterProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingEdgeThresholdParameterProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingNodeColorParameterProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingNodeFilteringProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingNodeRadiusParameterProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingNodeThresholdParameterProperty); MITK_REGISTER_ENUM_SUB_SERIALIZER(ConnectomicsRenderingSchemeProperty); #endif diff --git a/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp b/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp index da7a692..0611611 100644 --- a/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp +++ b/Modules/DiffusionCmdApps/FiberProcessing/FitFibersToImage.cpp @@ -1,321 +1,321 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include typedef itk::Point PointType4; typedef itk::Image< float, 4 > PeakImgType; /*! \brief Fits the tractogram to the input peak image by assigning a weight to each fiber (similar to https://doi.org/10.1016/j.neuroimage.2015.06.092). */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Fit Fibers To Image"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Assigns a weight to each fiber in order to optimally explain the input peak image"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i1", mitkDiffusionCommandLineParser::StringList, "Input tractograms:", "input tractograms (files or folder)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "i2", mitkDiffusionCommandLineParser::String, "Input image:", "input image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output root", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("max_iter", "", mitkDiffusionCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20); parser.addArgument("bundle_based", "", mitkDiffusionCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false); parser.addArgument("min_g", "", mitkDiffusionCommandLineParser::Float, "Min. g:", "lower termination threshold for gradient magnitude", 1e-5); parser.addArgument("lambda", "", mitkDiffusionCommandLineParser::Float, "Lambda:", "modifier for regularization", 1.0); parser.addArgument("save_res", "", mitkDiffusionCommandLineParser::Bool, "Save Residuals:", "save residual images", false); parser.addArgument("save_weights", "", mitkDiffusionCommandLineParser::Bool, "Save Weights:", "save fiber weights in a separate text file", false); parser.addArgument("filter_zero", "", mitkDiffusionCommandLineParser::Bool, "Filter Zero Weights:", "filter fibers with zero weight", false); parser.addArgument("filter_outliers", "", mitkDiffusionCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false); parser.addArgument("join_tracts", "", mitkDiffusionCommandLineParser::Bool, "Join output tracts:", "outout tracts are merged into a single tractogram", false); parser.addArgument("regu", "", mitkDiffusionCommandLineParser::String, "Regularization:", "MSM; Variance; VoxelVariance; Lasso; GroupLasso; GroupVariance; NONE", std::string("VoxelVariance")); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType fib_files = us::any_cast(parsedArgs["i1"]); std::string input_image_name = us::any_cast(parsedArgs["i2"]); std::string outRoot = us::any_cast(parsedArgs["o"]); bool single_fib = true; if (parsedArgs.count("bundle_based")) single_fib = !us::any_cast(parsedArgs["bundle_based"]); bool save_residuals = false; if (parsedArgs.count("save_res")) save_residuals = us::any_cast(parsedArgs["save_res"]); bool filter_zero = false; if (parsedArgs.count("filter_zero")) filter_zero = us::any_cast(parsedArgs["filter_zero"]); bool save_weights = false; if (parsedArgs.count("save_weights")) save_weights = us::any_cast(parsedArgs["save_weights"]); std::string regu = "VoxelVariance"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool join_tracts = false; if (parsedArgs.count("join_tracts")) join_tracts = us::any_cast(parsedArgs["join_tracts"]); int max_iter = 20; if (parsedArgs.count("max_iter")) max_iter = us::any_cast(parsedArgs["max_iter"]); float g_tol = 1e-5; if (parsedArgs.count("min_g")) g_tol = us::any_cast(parsedArgs["min_g"]); float lambda = 1.0; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); try { MITK_INFO << "Loading data"; mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, std::vector()); std::vector< std::string > fib_names; auto input_tracts = mitk::DiffusionDataIOHelper::load_fibs(fib_files, &fib_names); itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); mitk::BaseData::Pointer inputData = mitk::IOUtil::Load(input_image_name, &functor)[0].GetPointer(); mitk::Image::Pointer mitk_image = dynamic_cast(inputData.GetPointer()); mitk::PeakImage::Pointer mitk_peak_image = dynamic_cast(inputData.GetPointer()); if (mitk_peak_image.IsNotNull()) { typedef mitk::ImageToItk< mitk::PeakImage::ItkPeakImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitk_peak_image); caster->Update(); mitk::PeakImage::ItkPeakImageType::Pointer peak_image = caster->GetOutput(); fitter->SetPeakImage(peak_image); } else if (mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { fitter->SetDiffImage(mitk::DiffusionPropertyHelper::GetItkVectorImage(mitk_image)); mitk::TensorModel<>* model = new mitk::TensorModel<>(); model->SetBvalue(1000); model->SetDiffusivity1(0.0010); model->SetDiffusivity2(0.00015); model->SetDiffusivity3(0.00015); model->SetGradientList(mitk::DiffusionPropertyHelper::GetGradientContainer(mitk_image)); fitter->SetSignalModel(model); } else if (mitk_image->GetDimension()==3) { itk::FitFibersToImageFilter::DoubleImgType::Pointer scalar_image = itk::FitFibersToImageFilter::DoubleImgType::New(); mitk::CastToItkImage(mitk_image, scalar_image); fitter->SetScalarImage(scalar_image); } else { MITK_INFO << "Input image invalid. Valid options are peak image, 3D scalar image or raw diffusion-weighted image."; return EXIT_FAILURE; } fitter->SetTractograms(input_tracts); fitter->SetFitIndividualFibers(single_fib); fitter->SetMaxIterations(max_iter); fitter->SetGradientTolerance(g_tol); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); mitk::LocaleSwitch localeSwitch("C"); if (save_residuals && mitk_peak_image.IsNotNull()) { itk::ImageFileWriter< PeakImgType >::Pointer writer = itk::ImageFileWriter< PeakImgType >::New(); writer->SetInput(fitter->GetFittedImage()); writer->SetFileName(outRoot + "_fitted.nii.gz"); writer->Update(); writer->SetInput(fitter->GetResidualImage()); writer->SetFileName(outRoot + "_residual.nii.gz"); writer->Update(); writer->SetInput(fitter->GetOverexplainedImage()); writer->SetFileName(outRoot + "_overexplained.nii.gz"); writer->Update(); writer->SetInput(fitter->GetUnderexplainedImage()); writer->SetFileName(outRoot + "_underexplained.nii.gz"); writer->Update(); } else if (save_residuals && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(mitk_image)) { { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetFittedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_fitted_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetResidualImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_residual_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetOverexplainedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_overexplained_image.nii.gz"); } { mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( fitter->GetUnderexplainedImageDiff().GetPointer() ); mitk::DiffusionPropertyHelper::CopyProperties(mitk_image, outImage, true); mitk::DiffusionPropertyHelper::InitializeImage( outImage ); mitk::IOUtil::Save(outImage, "application/vnd.mitk.nii.gz", outRoot + "_underexplained_image.nii.gz"); } } else if (save_residuals) { itk::ImageFileWriter< itk::FitFibersToImageFilter::DoubleImgType >::Pointer writer = itk::ImageFileWriter< itk::FitFibersToImageFilter::DoubleImgType >::New(); writer->SetInput(fitter->GetFittedImageScalar()); writer->SetFileName(outRoot + "_fitted_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetResidualImageScalar()); writer->SetFileName(outRoot + "_residual_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetOverexplainedImageScalar()); writer->SetFileName(outRoot + "_overexplained_image.nii.gz"); writer->Update(); writer->SetInput(fitter->GetUnderexplainedImageScalar()); writer->SetFileName(outRoot + "_underexplained_image.nii.gz"); writer->Update(); } std::vector< mitk::FiberBundle::Pointer > output_tracts = fitter->GetTractograms(); if (!join_tracts) { for (unsigned int bundle=0; bundleFilterByWeights(0.0); if (fib->GetNumFibers()>0) { fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(fib, outRoot + name + "_fitted.fib"); if (save_weights) { - ofstream logfile; + std::ofstream logfile; logfile.open (outRoot + name + "_weights.txt"); for (unsigned int f=0; fGetNumFibers(); ++f) logfile << output_tracts.at(bundle)->GetFiberWeight(f) << "\n"; logfile.close(); } } else MITK_INFO << "Output contains no fibers!"; } } else { mitk::FiberBundle::Pointer out = mitk::FiberBundle::New(); out = out->AddBundles(output_tracts); if (filter_zero) out = out->FilterByWeights(0.0); if (out->GetNumFibers()>0) { out->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out, outRoot + "_fitted.fib"); if (save_weights) { - ofstream logfile; + std::ofstream logfile; logfile.open (outRoot + "_weights.txt"); for (unsigned int f=0; fGetNumFibers(); ++f) logfile << out->GetFiberWeight(f) << "\n"; logfile.close(); } } else MITK_INFO << "Output contains no fibers!"; } } catch (const itk::ExceptionObject& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/FiberProcessing/PrintFiberStatistics.cpp b/Modules/DiffusionCmdApps/FiberProcessing/PrintFiberStatistics.cpp index a166121..1a716d6 100644 --- a/Modules/DiffusionCmdApps/FiberProcessing/PrintFiberStatistics.cpp +++ b/Modules/DiffusionCmdApps/FiberProcessing/PrintFiberStatistics.cpp @@ -1,101 +1,101 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Join multiple tractograms */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Fiber Statistics"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Output statistics about fiber lengths, weights, etc. to a file."); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input:", "input tractograms", us::Any(), false); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output File:", "output file", us::Any(), false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType inFibs = us::any_cast(parsedArgs["i"]); std::string outfile = us::any_cast(parsedArgs["o"]); try { std::vector< std::string > fib_names; auto input_tracts = mitk::DiffusionDataIOHelper::load_fibs(inFibs, &fib_names); - ofstream statistics_file; + std::ofstream statistics_file; statistics_file.open (outfile); int c=0; for (auto tract : input_tracts) { MITK_INFO << "Writing statistics of bunlde " << fib_names.at(c); statistics_file << "File: " << fib_names.at(c) << std::endl; tract->PrintSelf(statistics_file, 2); statistics_file << "-------------------------------------------------------------------------\n\n" << std::endl; ++c; } statistics_file.close(); } catch (const itk::ExceptionObject& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp b/Modules/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp index d64e8e9..b227ea6 100644 --- a/Modules/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp +++ b/Modules/DiffusionCmdApps/TractographyEvaluation/AnchorConstrainedPlausibility.cpp @@ -1,447 +1,447 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include typedef itk::Point PointType4; typedef mitk::PeakImage::ItkPeakImageType PeakImgType; typedef itk::Image< unsigned char, 3 > ItkUcharImageType; /*! \brief Score input candidate tracts using ACP analysis */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Anchor Constrained Plausibility"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription("Score input candidate tracts using ACP analysis"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "a", mitkDiffusionCommandLineParser::String, "Anchor tractogram:", "anchor tracts in one tractogram file", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "p", mitkDiffusionCommandLineParser::String, "Input peaks:", "input peak image", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "c", mitkDiffusionCommandLineParser::StringList, "Candidates:", "Folder(s) or file list of candidate tracts", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output folder:", "output folder", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("reference_mask_folders", "", mitkDiffusionCommandLineParser::StringList, "Reference Mask Folder(s):", "Folder(s) or file list containing reference tract masks for accuracy evaluation", true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("reference_peaks_folders", "", mitkDiffusionCommandLineParser::StringList, "Reference Peaks Folder(s):", "Folder(s) or file list containing reference peak images for accuracy evaluation", true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("mask", "", mitkDiffusionCommandLineParser::String, "Mask image:", "scoring is only performed inside the mask image", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("greedy_add", "", mitkDiffusionCommandLineParser::Bool, "Greedy:", "if enabled, the candidate tracts are not jointly fitted to the residual image but one after the other employing a greedy scheme", false); parser.addArgument("lambda", "", mitkDiffusionCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1); parser.addArgument("filter_outliers", "", mitkDiffusionCommandLineParser::Bool, "Filter outliers:", "perform second optimization run with an upper weight bound based on the first weight estimation (99% quantile)", false); parser.addArgument("regu", "", mitkDiffusionCommandLineParser::String, "Regularization:", "MSM; Variance; VoxelVariance; Lasso; GroupLasso; GroupVariance; NONE", std::string("NONE")); parser.addArgument("use_num_streamlines", "", mitkDiffusionCommandLineParser::Bool, "Use number of streamlines as score:", "Don't fit candidates, simply use number of streamlines per candidate as score", false); parser.addArgument("use_weights", "", mitkDiffusionCommandLineParser::Bool, "Use input weights as score:", "Don't fit candidates, simply use first input streamline weight per candidate as score", false); parser.addArgument("filter_zero_weights", "", mitkDiffusionCommandLineParser::Bool, "Filter zero-weights", "Remove streamlines with weight 0 from candidates", false); parser.addArgument("flipx", "", mitkDiffusionCommandLineParser::Bool, "Flip x", "flip along x-axis", false); parser.addArgument("flipy", "", mitkDiffusionCommandLineParser::Bool, "Flip y", "flip along y-axis", false); parser.addArgument("flipz", "", mitkDiffusionCommandLineParser::Bool, "Flip z", "flip along z-axis", false); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string peak_file_name = us::any_cast(parsedArgs["p"]); std::string out_folder = us::any_cast(parsedArgs["o"]); mitkDiffusionCommandLineParser::StringContainerType candidate_tract_folders = us::any_cast(parsedArgs["c"]); if (!out_folder.empty() && out_folder.back() != '/') out_folder += "/"; bool greedy_add = false; if (parsedArgs.count("greedy_add")) greedy_add = us::any_cast(parsedArgs["greedy_add"]); float lambda = 0.1f; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); bool filter_outliers = false; if (parsedArgs.count("filter_outliers")) filter_outliers = us::any_cast(parsedArgs["filter_outliers"]); bool filter_zero_weights = false; if (parsedArgs.count("filter_zero_weights")) filter_zero_weights = us::any_cast(parsedArgs["filter_zero_weights"]); std::string mask_file = ""; if (parsedArgs.count("mask")) mask_file = us::any_cast(parsedArgs["mask"]); mitkDiffusionCommandLineParser::StringContainerType reference_mask_files_folders; if (parsedArgs.count("reference_mask_folders")) reference_mask_files_folders = us::any_cast(parsedArgs["reference_mask_folders"]); mitkDiffusionCommandLineParser::StringContainerType reference_peaks_files_folders; if (parsedArgs.count("reference_peaks_folders")) reference_peaks_files_folders = us::any_cast(parsedArgs["reference_peaks_folders"]); std::string regu = "NONE"; if (parsedArgs.count("regu")) regu = us::any_cast(parsedArgs["regu"]); bool use_weights = false; if (parsedArgs.count("use_weights")) use_weights = us::any_cast(parsedArgs["use_weights"]); bool use_num_streamlines = false; if (parsedArgs.count("use_num_streamlines")) use_num_streamlines = us::any_cast(parsedArgs["use_num_streamlines"]); bool flipx = false; if (parsedArgs.count("flipx")) flipx = us::any_cast(parsedArgs["flipx"]); bool flipy = false; if (parsedArgs.count("flipy")) flipy = us::any_cast(parsedArgs["flipy"]); bool flipz = false; if (parsedArgs.count("flipz")) flipz = us::any_cast(parsedArgs["flipz"]); try { itk::TimeProbe clock; clock.Start(); if (!ist::PathExists(out_folder)) { MITK_INFO << "Creating output directory"; ist::MakeDirectory(out_folder); } MITK_INFO << "Loading data"; // Load mask file. Fit is only performed inside the mask MITK_INFO << "Loading mask image"; auto mask = mitk::DiffusionDataIOHelper::load_itk_image(mask_file); // Load masks covering the true positives for evaluation purposes MITK_INFO << "Loading reference peaks and masks"; std::vector< std::string > anchor_mask_files; auto reference_masks = mitk::DiffusionDataIOHelper::load_itk_images(reference_mask_files_folders, &anchor_mask_files); auto reference_peaks = mitk::DiffusionDataIOHelper::load_itk_images(reference_peaks_files_folders); // Load peak image MITK_INFO << "Loading peak image"; auto peak_image = mitk::DiffusionDataIOHelper::load_itk_image(peak_file_name); // Load all candidate tracts MITK_INFO << "Loading candidate tracts"; std::vector< std::string > candidate_tract_files; auto input_candidates = mitk::DiffusionDataIOHelper::load_fibs(candidate_tract_folders, &candidate_tract_files); if (flipx || flipy || flipz) { itk::FlipPeaksFilter< float >::Pointer flipper = itk::FlipPeaksFilter< float >::New(); flipper->SetInput(peak_image); flipper->SetFlipX(flipx); flipper->SetFlipY(flipy); flipper->SetFlipZ(flipz); flipper->Update(); peak_image = flipper->GetOutput(); } mitk::LocaleSwitch localeSwitch("C"); itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New(); - ofstream logfile; + std::ofstream logfile; logfile.open (out_folder + "scores.txt"); double rmse = 0.0; int iteration = 0; std::string name = "NOANCHOR"; if (parsedArgs.count("a")) { // Load reference tractogram consisting of all known tracts std::string anchors_file = us::any_cast(parsedArgs["a"]); mitk::FiberBundle::Pointer anchor_tractogram = mitk::IOUtil::Load(anchors_file); if ( !(anchor_tractogram.IsNull() || anchor_tractogram->GetNumFibers()==0) ) { // Fit known tracts to peak image to obtain underexplained image MITK_INFO << "Fit anchor tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetTractograms({anchor_tractogram}); fitter->SetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetPeakImage(peak_image); fitter->SetVerbose(true); fitter->SetMaskImage(mask); fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); rmse = fitter->GetRMSE(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); name = ist::GetFilenameWithoutExtension(anchors_file); mitk::FiberBundle::Pointer anchor_tracts = fitter->GetTractograms().at(0); anchor_tracts->SetFiberColors(255,255,255); mitk::IOUtil::Save(anchor_tracts, out_folder + boost::lexical_cast(static_cast(100000*rms_diff[0])) + "_" + name + ".fib"); logfile << name << " " << setprecision(5) << rms_diff[0] << "\n"; peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_" + name + ".nii.gz"); peak_image_writer->Update(); } } if (use_weights || use_num_streamlines) { MITK_INFO << "Using tract weights as scores"; unsigned int c = 0; for (auto fib : input_candidates) { int mod = 1; float score = 0; if (use_weights) { score = fib->GetFiberWeight(0); mod = 100000; } else if (use_num_streamlines) score = fib->GetNumFibers(); fib->ColorFibersByOrientation(); std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast(static_cast(mod*score)) + "_" + bundle_name + ".fib"); unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } float weight_sum = 0; for (unsigned int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << bundle_name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } } else if (!greedy_add) { MITK_INFO << "Fit candidate tracts"; itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); fitter->SetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(true); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); fitter->SetTractograms(input_candidates); fitter->SetFitIndividualFibers(true); if (regu=="MSM") fitter->SetRegularization(VnlCostFunction::REGU::MSM); else if (regu=="Variance") fitter->SetRegularization(VnlCostFunction::REGU::VARIANCE); else if (regu=="Lasso") fitter->SetRegularization(VnlCostFunction::REGU::LASSO); else if (regu=="VoxelVariance") fitter->SetRegularization(VnlCostFunction::REGU::VOXEL_VARIANCE); else if (regu=="GroupLasso") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_LASSO); else if (regu=="GroupVariance") fitter->SetRegularization(VnlCostFunction::REGU::GROUP_VARIANCE); else if (regu=="NONE") fitter->SetRegularization(VnlCostFunction::REGU::NONE); fitter->Update(); vnl_vector rms_diff = fitter->GetRmsDiffPerBundle(); unsigned int c = 0; for (auto fib : input_candidates) { std::string bundle_name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(c)); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect if (filter_zero_weights) fib = fib->FilterByWeights(0); mitk::IOUtil::Save(fib, out_folder + boost::lexical_cast((int)(100000*rms_diff[c])) + "_" + bundle_name + ".fib"); unsigned int num_voxels = 0; { itk::TractDensityImageFilter< ItkUcharImageType >::Pointer masks_filter = itk::TractDensityImageFilter< ItkUcharImageType >::New(); masks_filter->SetInputImage(mask); masks_filter->SetBinaryOutput(true); masks_filter->SetFiberBundle(fib); masks_filter->SetUseImageGeometry(true); masks_filter->Update(); num_voxels = masks_filter->GetNumCoveredVoxels(); } float weight_sum = 0; for (unsigned int i=0; iGetNumFibers(); i++) weight_sum += fib->GetFiberWeight(i); std::cout.rdbuf (old); // <-- restore logfile << bundle_name << " " << setprecision(5) << rms_diff[c] << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; ++c; } mitk::FiberBundle::Pointer out_fib = mitk::FiberBundle::New(); out_fib = out_fib->AddBundles(input_candidates); out_fib->ColorFibersByFiberWeights(false, true); mitk::IOUtil::Save(out_fib, out_folder + "AllCandidates.fib"); peak_image = fitter->GetUnderexplainedImage(); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + "Residual_AllCandidates.nii.gz"); peak_image_writer->Update(); } else { MITK_INFO << "RMSE: " << setprecision(5) << rmse; // fitter->SetPeakImage(peak_image); // Iteratively add candidate bundles in a greedy manner while (!input_candidates.empty()) { double next_rmse = rmse; mitk::FiberBundle::Pointer best_candidate = nullptr; PeakImgType::Pointer best_candidate_peak_image = nullptr; for (unsigned int i=0; iSetLambda(static_cast(lambda)); fitter->SetFilterOutliers(filter_outliers); fitter->SetVerbose(false); fitter->SetPeakImage(peak_image); fitter->SetMaskImage(mask); // ****************************** fitter->SetTractograms({input_candidates.at(i)}); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect fitter->Update(); std::cout.rdbuf (old); // <-- restore double candidate_rmse = fitter->GetRMSE(); if (candidate_rmseGetTractograms().at(0); best_candidate_peak_image = fitter->GetUnderexplainedImage(); } } if (best_candidate.IsNull()) break; // fitter->SetPeakImage(peak_image); peak_image = best_candidate_peak_image; unsigned int i=0; std::vector< mitk::FiberBundle::Pointer > remaining_candidates; std::vector< std::string > remaining_candidate_files; for (auto fib : input_candidates) { if (fib!=best_candidate) { remaining_candidates.push_back(fib); remaining_candidate_files.push_back(candidate_tract_files.at(i)); } else name = ist::GetFilenameWithoutExtension(candidate_tract_files.at(i)); ++i; } input_candidates = remaining_candidates; candidate_tract_files = remaining_candidate_files; iteration++; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect // Save winning candidate if (filter_zero_weights) best_candidate = best_candidate->FilterByWeights(0); mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast(iteration) + "_" + name + ".fib"); peak_image_writer->SetInput(peak_image); peak_image_writer->SetFileName(out_folder + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); peak_image_writer->Update(); std::cout.rdbuf (old); // <-- restore // logfile << name << " " << setprecision(5) << score << " " << num_voxels << " " << fib->GetNumFibers() << " " << weight_sum << "\n"; } } clock.Stop(); int h = static_cast(clock.GetTotal()/3600); int m = (static_cast(clock.GetTotal())%3600)/60; int s = static_cast(clock.GetTotal())%60; MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s"; logfile.close(); } catch (const itk::ExceptionObject& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp b/Modules/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp index e106074..67879d3 100644 --- a/Modules/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp +++ b/Modules/DiffusionCmdApps/TractographyEvaluation/GetOverlappingTracts.cpp @@ -1,167 +1,167 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Get Overlapping Tracts"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setDescription("Find tracts that overlap with the reference masks or tracts"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input:", "input tractograms (.fib/.trk/.tck/.dcm)", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output Folder:", "move input tracts that do/don't overlap here", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("reference", "", mitkDiffusionCommandLineParser::StringList, "Reference:", "reference tractograms or mask images", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("overlap_fraction", "", mitkDiffusionCommandLineParser::Float, "Overlap fraction:", "", 0.9); parser.addArgument("use_any_overlap", "", mitkDiffusionCommandLineParser::Bool, "Use any overlap:", "Don't find maximum overlap but use first overlap larger threshold"); parser.addArgument("dont_save_tracts", "", mitkDiffusionCommandLineParser::Bool, "Don't save tracts:", "if true, only text files documenting the overlaps are saved and no tract files are copied"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType input = us::any_cast(parsedArgs["i"]); mitkDiffusionCommandLineParser::StringContainerType reference = us::any_cast(parsedArgs["reference"]); std::string out_folder = us::any_cast(parsedArgs["o"]); bool use_any_overlap = false; if (parsedArgs.count("use_any_overlap")) use_any_overlap = us::any_cast(parsedArgs["use_any_overlap"]); bool dont_save_tracts = false; if (parsedArgs.count("dont_save_tracts")) dont_save_tracts = us::any_cast(parsedArgs["dont_save_tracts"]); float overlap_threshold = 0.9; if (parsedArgs.count("overlap_fraction")) overlap_threshold = us::any_cast(parsedArgs["overlap_fraction"]); try { MITK_INFO << "Loading references"; std::vector< std::string > reference_names; auto masks = mitk::DiffusionDataIOHelper::load_itk_images(reference, &reference_names); auto reference_fibs = mitk::DiffusionDataIOHelper::load_fibs(reference, &reference_names); std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect itk::TractDensityImageFilter< ItkFloatImgType >::Pointer filter = itk::TractDensityImageFilter< ItkFloatImgType >::New(); filter->SetUpsamplingFactor(0.25); filter->SetBinaryOutput(true); for (auto fib : reference_fibs) { filter->SetFiberBundle(fib); filter->Update(); masks.push_back(filter->GetOutput()); } std::cout.rdbuf (old); // <-- restore MITK_INFO << "Loading input tractograms"; std::vector< std::string > input_names; auto input_fibs = mitk::DiffusionDataIOHelper::load_fibs(input, &input_names); MITK_INFO << "Finding overlaps"; - ofstream logfile; + std::ofstream logfile; logfile.open (out_folder + "Overlaps.txt"); - ofstream logfile2; + std::ofstream logfile2; logfile2.open (out_folder + "AllOverlaps.txt"); boost::progress_display disp(input.size()); unsigned int c = 0; for (auto fib : input_fibs) { ++disp; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect bool is_overlapping = false; float overlap = 0; float max_overlap = 0; std::string max_ref = "-"; int i = 0; std::string overlap_string = ist::GetFilenameWithoutExtension(input_names.at(c)); for (auto m : masks) { overlap = fib->GetOverlap(m); if (overlap>max_overlap) { max_overlap = overlap; max_ref = ist::GetFilenameWithoutExtension(reference_names.at(i)); } if (use_any_overlap && overlap>=overlap_threshold) break; overlap_string += " " + ist::GetFilenameWithoutExtension(reference_names.at(i)) + " " + boost::lexical_cast(overlap); ++i; } if (overlap>=overlap_threshold) is_overlapping = true; logfile << ist::GetFilenameWithoutExtension(input_names.at(c)) << " - " << max_ref << ": " << boost::lexical_cast(max_overlap) << "\n"; logfile2 << overlap_string << "\n"; if (!dont_save_tracts && is_overlapping) ist::CopyAFile(input_names.at(c), out_folder + ist::GetFilenameName(input_names.at(c))); std::cout.rdbuf (old); // <-- restore ++c; } logfile.close(); logfile2.close(); } catch (const itk::ExceptionObject& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp b/Modules/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp index b9d4b62..806f899 100644 --- a/Modules/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp +++ b/Modules/DiffusionCmdApps/TractographyEvaluation/PeaksAngularError.cpp @@ -1,190 +1,190 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include "mitkDiffusionCommandLineParser.h" #include #include #include #include #include #include #include #include #include #include typedef itk::Image< unsigned char, 3 > ItkUcharImageType; /*! \brief Calculate angular error between two sets of directions stored in multiple 3D vector images where each pixel corresponds to a vector (itk::Image< itk::Vector< float, 3>, 3 >) */ int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setArgumentPrefix("--", "-"); parser.addArgument("test", "", mitkDiffusionCommandLineParser::StringList, "Test images", "test direction images", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("reference", "", mitkDiffusionCommandLineParser::StringList, "Reference images", "reference direction images", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output folder", "output folder", us::Any(), false, false, false, mitkDiffusionCommandLineParser::Output); parser.addArgument("masks", "", mitkDiffusionCommandLineParser::StringList, "Mask(s)", "mask image(s)", us::Any(), true, false, false, mitkDiffusionCommandLineParser::Input); parser.addArgument("verbose", "", mitkDiffusionCommandLineParser::Bool, "Verbose", "output error images"); parser.addArgument("ignore_test", "", mitkDiffusionCommandLineParser::Bool, "Ignore missing test", "don't increase error if no test directions are found"); parser.addArgument("ignore_ref", "", mitkDiffusionCommandLineParser::Bool, "Ignore ignore missing ref", "don't increase error if no ref directions are found"); parser.setCategory("Fiber Tracking Evaluation"); parser.setTitle("Peaks Angular Error"); parser.setDescription("Calculate angular error between two sets of peak images (1-1 correspondence)"); parser.setContributor("MIC"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType testImages = us::any_cast(parsedArgs["test"]); mitkDiffusionCommandLineParser::StringContainerType referenceImages = us::any_cast(parsedArgs["reference"]); mitkDiffusionCommandLineParser::StringContainerType maskImages; if (parsedArgs.count("masks")) maskImages = us::any_cast(parsedArgs["masks"]); std::string outRoot = us::any_cast(parsedArgs["o"]); bool verbose = false; if (parsedArgs.count("verbose")) verbose = us::any_cast(parsedArgs["verbose"]); bool ignore_test = false; if (parsedArgs.count("ignore_test")) ignore_test = us::any_cast(parsedArgs["ignore_test"]); bool ignore_ref = false; if (parsedArgs.count("ignore_ref")) ignore_ref = us::any_cast(parsedArgs["ignore_ref"]); try { typedef itk::ComparePeakImagesFilter< float > EvaluationFilterType; std::vector test_names; auto test_images = mitk::DiffusionDataIOHelper::load_itk_images(testImages, &test_names); // load reference directions std::vector ref_names; auto ref_images = mitk::DiffusionDataIOHelper::load_itk_images(referenceImages, &ref_names); // load/create mask image auto itkMaskImages = mitk::DiffusionDataIOHelper::load_itk_images(maskImages); if (test_images.size()!=ref_images.size()) mitkThrow() << "Matching number of test and reference image required!"; for (unsigned int i=0; iSetTestImage(test_images.at(i)); evaluationFilter->SetReferenceImage(ref_images.at(i)); if (iSetMaskImage(itkMaskImages.at(i)); evaluationFilter->SetIgnoreMissingTestDirections(ignore_test); evaluationFilter->SetIgnoreMissingRefDirections(ignore_ref); evaluationFilter->Update(); std::string ref_name = ist::GetFilenameWithoutExtension(ref_names.at(i)); std::string test_name = ist::GetFilenameWithoutExtension(test_names.at(i)); if (verbose) { mitk::LocaleSwitch localeSwitch("C"); EvaluationFilterType::OutputImageType::Pointer angularErrorImage = evaluationFilter->GetOutput(0); EvaluationFilterType::OutputImageType::Pointer lengthErrorImage = evaluationFilter->GetOutput(1); typedef itk::ImageFileWriter< EvaluationFilterType::OutputImageType > WriterType; { WriterType::Pointer writer = WriterType::New(); std::string outfilename = outRoot; outfilename.append(ref_name + "_" + test_name + "_AngularError.nii.gz"); writer->SetFileName(outfilename.c_str()); writer->SetInput(angularErrorImage); writer->Update(); } { WriterType::Pointer writer = WriterType::New(); std::string outfilename = outRoot; outfilename.append(ref_name + "_" + test_name + "_LengthError.nii.gz"); writer->SetFileName(outfilename.c_str()); writer->SetInput(lengthErrorImage); writer->Update(); } } std::string logFile = outRoot; logFile.append("AngularErrors.csv"); bool add_header = true; if (ist::FileExists(logFile, true)) add_header = false; - ofstream file; + std::ofstream file; file.open (logFile.c_str(), std::fstream::app); std::string sens; if (add_header) sens.append("Test,Reference,Mean,Median,Maximum,Minimum,Stdev\n"); sens.append(test_name); sens.append(","); sens.append(ref_name); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMeanAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMedianAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMaxAngularError())); sens.append(","); sens.append(boost::lexical_cast(evaluationFilter->GetMinAngularError())); sens.append(","); sens.append(boost::lexical_cast(std::sqrt(evaluationFilter->GetVarAngularError()))); sens.append("\n"); std::cout << sens; file << sens; file.close(); } } catch (const itk::ExceptionObject& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (std::exception& e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; } diff --git a/Modules/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp b/Modules/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp index 54666d3..dd32fc8 100644 --- a/Modules/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp +++ b/Modules/DiffusionCmdApps/TractographyEvaluation/ReferenceSimilarity.cpp @@ -1,149 +1,149 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include typedef itk::Image< unsigned char, 3 > ItkUcharImageType; typedef itk::Image< float, 4 > ItkPeakImgType; int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Reference Similarity"); parser.setCategory("Fiber Tracking Evaluation"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkDiffusionCommandLineParser::StringList, "Input Tracts:", "input tracts folder", us::Any(), false); parser.addArgument("reference_tracts", "", mitkDiffusionCommandLineParser::StringList, "", "", us::Any(), false); parser.addArgument("reference_masks", "", mitkDiffusionCommandLineParser::StringList, "", "", us::Any(), false); parser.addArgument("reference_peaks", "", mitkDiffusionCommandLineParser::StringList, "", "", us::Any(), false); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "", "", us::Any(), false); parser.addArgument("fiber_points", "", mitkDiffusionCommandLineParser::Int, "Fiber points:", "", 20); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string out_folder = us::any_cast(parsedArgs["o"]); mitkDiffusionCommandLineParser::StringContainerType input_tract_files = us::any_cast(parsedArgs["i"]); mitkDiffusionCommandLineParser::StringContainerType reference_tract_files = us::any_cast(parsedArgs["reference_tracts"]); mitkDiffusionCommandLineParser::StringContainerType reference_mask_files = us::any_cast(parsedArgs["reference_masks"]); mitkDiffusionCommandLineParser::StringContainerType reference_peak_files = us::any_cast(parsedArgs["reference_peaks"]); int fiber_points = 20; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); try { std::vector input_tract_names; std::vector ref_tract_names; std::vector< mitk::FiberBundle::Pointer > input_tracts = mitk::DiffusionDataIOHelper::load_fibs(input_tract_files, &input_tract_names); std::vector< mitk::FiberBundle::Pointer > reference_tracts = mitk::DiffusionDataIOHelper::load_fibs(reference_tract_files, &ref_tract_names); std::vector< ItkUcharImageType::Pointer > reference_masks = mitk::DiffusionDataIOHelper::load_itk_images(reference_mask_files); std::vector< ItkPeakImgType::Pointer > reference_peaks = mitk::DiffusionDataIOHelper::load_itk_images(reference_peak_files); MITK_INFO << "Calculating distances"; itk::TractDistanceFilter::Pointer distance_calculator = itk::TractDistanceFilter::New(); distance_calculator->SetNumPoints(fiber_points); distance_calculator->SetTracts1(input_tracts); distance_calculator->SetTracts2(reference_tracts); distance_calculator->SetMetrics({new mitk::ClusteringMetricEuclideanMean()}); distance_calculator->Update(); auto distances = distance_calculator->GetAllDistances(); vnl_matrix voxel_overlap; voxel_overlap.set_size(input_tracts.size(), reference_tracts.size()); vnl_matrix dir_overlap; dir_overlap.set_size(input_tracts.size(), reference_tracts.size()); MITK_INFO << "Calculating overlap"; boost::progress_display disp(input_tracts.size()*reference_tracts.size()); int r=0; for (auto fib : input_tracts) { int c=0; for (auto ref_mask : reference_masks) { // ++disp; // std::streambuf *old = cout.rdbuf(); // <-- save // std::stringstream ss; // std::cout.rdbuf (ss.rdbuf()); // <-- redirect float overlap = 0; float directional_overlap = 0; std::tie(directional_overlap, overlap) = fib->GetDirectionalOverlap(ref_mask, reference_peaks.at(c)); voxel_overlap[r][c] = overlap; dir_overlap[r][c] = directional_overlap; // std::cout.rdbuf (old); // <-- restore ++c; } ++r; } - ofstream logfile; + std::ofstream logfile; logfile.open(out_folder + "ref_tract_names.txt"); for (unsigned int i=0; i #include #include #include #include #include #include #include #include int main(int argc, char* argv[]) { mitkDiffusionCommandLineParser parser; parser.setTitle("Tract Distance"); parser.setCategory("Fiber Processing"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i1", mitkDiffusionCommandLineParser::StringList, "Input tracts 1:", "input tracts 1", us::Any(), false); parser.addArgument("", "i2", mitkDiffusionCommandLineParser::StringList, "Input tracts 2:", "input tracts 2", us::Any(), false); parser.addArgument("", "o", mitkDiffusionCommandLineParser::String, "Output:", "output logfile", us::Any(), false); parser.addArgument("fiber_points", "", mitkDiffusionCommandLineParser::Int, "Fiber points:", "", 12); parser.addArgument("metrics", "", mitkDiffusionCommandLineParser::StringList, "Metrics:", "EU_MEAN (default), EU_STD, EU_MAX"); parser.addArgument("metric_weights", "", mitkDiffusionCommandLineParser::StringList, "Metric weights:", "add one float weight for each used metric"); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkDiffusionCommandLineParser::StringContainerType t1_folder = us::any_cast(parsedArgs["i1"]); mitkDiffusionCommandLineParser::StringContainerType t2_folder = us::any_cast(parsedArgs["i2"]); std::string out_file = us::any_cast(parsedArgs["o"]); int fiber_points = 12; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); std::vector< std::string > metric_strings = {"EU_MEAN"}; if (parsedArgs.count("metrics")) metric_strings = us::any_cast(parsedArgs["metrics"]); std::vector< std::string > metric_weights = {"1.0"}; if (parsedArgs.count("metric_weights")) metric_weights = us::any_cast(parsedArgs["metric_weights"]); if (metric_strings.size()!=metric_weights.size()) { MITK_INFO << "Each metric needs an associated metric weight!"; return EXIT_FAILURE; } try { std::vector t1_files; std::vector< mitk::FiberBundle::Pointer > tractograms1 = mitk::DiffusionDataIOHelper::load_fibs(t1_folder, &t1_files); std::vector t2_files; std::vector< mitk::FiberBundle::Pointer > tractograms2 = mitk::DiffusionDataIOHelper::load_fibs(t2_folder, &t2_files); MITK_INFO << "Loaded " << tractograms1.size() << " source tractograms."; MITK_INFO << "Loaded " << tractograms2.size() << " target tractograms."; itk::TractDistanceFilter::Pointer distance_calculator = itk::TractDistanceFilter::New(); distance_calculator->SetNumPoints(fiber_points); distance_calculator->SetTracts1(tractograms1); distance_calculator->SetTracts2(tractograms2); std::vector< mitk::ClusteringMetric* > metrics; int mc = 0; for (auto m : metric_strings) { float w = boost::lexical_cast(metric_weights.at(mc)); MITK_INFO << "Metric: " << m << " (w=" << w << ")"; if (m=="EU_MEAN") metrics.push_back({new mitk::ClusteringMetricEuclideanMean()}); else if (m=="EU_STD") metrics.push_back({new mitk::ClusteringMetricEuclideanStd()}); else if (m=="EU_MAX") metrics.push_back({new mitk::ClusteringMetricEuclideanMax()}); metrics.back()->SetScale(w); mc++; } if (metrics.empty()) { MITK_INFO << "No metric selected!"; return EXIT_FAILURE; } distance_calculator->SetMetrics(metrics); distance_calculator->Update(); MITK_INFO << "Distances:"; auto distances = distance_calculator->GetMinDistances(); auto indices = distance_calculator->GetMinIndices(); - ofstream logfile; + std::ofstream logfile; logfile.open (out_file); for (unsigned int i=0; i #include #include #include #include #include "mitkPixelType.h" #include "itkImageRegionIterator.h" #include "itkImageFileReader.h" #include #include #include template< class D, class T > mitk::TeemDiffusionTensor3DReconstructionImageFilter ::TeemDiffusionTensor3DReconstructionImageFilter(): m_EstimateErrorImage(false),m_Sigma(-19191919), m_EstimationMethod(TeemTensorEstimationMethodsLLS), m_NumIterations(1),m_ConfidenceThreshold(-19191919.0), m_ConfidenceFuzzyness(0.0),m_MinPlausibleValue(1.0) { } template< class D, class T > mitk::TeemDiffusionTensor3DReconstructionImageFilter ::~TeemDiffusionTensor3DReconstructionImageFilter() { } void file_replace(std::string filename, std::string what, std::string with) { - ofstream myfile2; + std::ofstream myfile2; std::locale C("C"); std::locale originalLocale2 = myfile2.getloc(); myfile2.imbue(C); char filename2[512]; sprintf(filename2, "%s2",filename.c_str()); myfile2.open (filename2); std::string line; - ifstream myfile (filename.c_str()); + std::ifstream myfile (filename.c_str()); std::locale originalLocale = myfile.getloc(); myfile.imbue(C); if (myfile.is_open()) { while (! myfile.eof() ) { getline (myfile,line); itksys::SystemTools::ReplaceString(line,what.c_str(),with.c_str()); myfile2 << line << std::endl; } myfile.close(); } myfile2.close(); itksys::SystemTools::RemoveFile(filename.c_str()); rename(filename2,filename.c_str()); myfile.imbue( originalLocale ); myfile2.imbue( originalLocale2 ); } // do the work template< class D, class T > void mitk::TeemDiffusionTensor3DReconstructionImageFilter ::Update() { itk::Statistics::MersenneTwisterRandomVariateGenerator::Pointer randgen = itk::Statistics::MersenneTwisterRandomVariateGenerator::New(); randgen->SetSeed(); // save input image to nrrd file in temp-folder char filename[512]; int random_integer = randgen->GetIntegerVariate(); sprintf( filename, "dwi_%d.nhdr",random_integer); try { mitk::IOUtil::Save(m_Input, filename); } catch (const itk::ExceptionObject& e) { std::cout << e.GetDescription() << std::endl; } file_replace(filename,"vector","list"); // build up correct command from input params char command[4096]; sprintf( command, "tend estim -i %s -B kvp -o tensors_%d.nhdr -knownB0 true", filename, random_integer); //m_DiffusionImages; if(m_EstimateErrorImage) { sprintf( command, "%s -ee error_image_%d.nhdr", command, random_integer); } if(m_Sigma != -19191919) { sprintf( command, "%s -sigma %f", command, m_Sigma); } switch(m_EstimationMethod) { case TeemTensorEstimationMethodsLLS: sprintf( command, "%s -est lls", command); break; case TeemTensorEstimationMethodsMLE: sprintf( command, "%s -est mle", command); break; case TeemTensorEstimationMethodsNLS: sprintf( command, "%s -est nls", command); break; case TeemTensorEstimationMethodsWLS: sprintf( command, "%s -est wls", command); break; } sprintf( command, "%s -wlsi %d", command, m_NumIterations); if(m_ConfidenceThreshold != -19191919.0) { sprintf( command, "%s -t %f", command, m_ConfidenceThreshold); } sprintf( command, "%s -soft %f", command, m_ConfidenceFuzzyness); sprintf( command, "%s -mv %f", command, m_MinPlausibleValue); // call tend estim command std::cout << "Calling <" << command << ">" << std::endl; int success = system(command); if(!success) { MITK_ERROR << "system command could not be called!"; } remove(filename); sprintf( filename, "dwi_%d.raw", random_integer); remove(filename); // change kind from tensor to vector sprintf( filename, "tensors_%d.nhdr", random_integer); file_replace(filename,"3D-masked-symmetric-matrix","vector"); mitk::LocaleSwitch localeSwitch("C"); // read result as mitk::Image and provide it in m_Output typedef itk::ImageFileReader FileReaderType; typename FileReaderType::Pointer reader = FileReaderType::New(); reader->SetFileName(filename); reader->Update(); typename VectorImageType::Pointer vecImage = reader->GetOutput(); remove(filename); sprintf( filename, "tensors_%d.raw", random_integer); remove(filename); typename ItkTensorImageType::Pointer itkTensorImage = ItkTensorImageType::New(); itkTensorImage->SetSpacing( vecImage->GetSpacing() ); // Set the image spacing itkTensorImage->SetOrigin( vecImage->GetOrigin() ); // Set the image origin itkTensorImage->SetDirection( vecImage->GetDirection() ); // Set the image direction itkTensorImage->SetLargestPossibleRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->SetBufferedRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->SetRequestedRegion( vecImage->GetLargestPossibleRegion() ); itkTensorImage->Allocate(); itk::ImageRegionIterator it(vecImage, vecImage->GetLargestPossibleRegion()); itk::ImageRegionIterator it2(itkTensorImage, itkTensorImage->GetLargestPossibleRegion()); it2 = it2.Begin(); //#pragma omp parallel private (it) { for(it=it.Begin();!it.IsAtEnd(); ++it, ++it2) { //#pragma omp single nowait { VectorType vec = it.Get(); TensorType tensor; for(int i=1;i<7;i++) tensor[i-1] = vec[i] * vec[0]; it2.Set( tensor ); } } // end for } // end ompparallel m_OutputItk = mitk::TensorImage::New(); m_OutputItk->InitializeByItk(itkTensorImage.GetPointer()); m_OutputItk->SetVolume( itkTensorImage->GetBufferPointer() ); // in case: read resulting error-image and provide it in m_ErrorImage if(m_EstimateErrorImage) { // open error image here } } diff --git a/Modules/DiffusionCore/CMakeLists.txt b/Modules/DiffusionCore/CMakeLists.txt index cdab554..ded455b 100644 --- a/Modules/DiffusionCore/CMakeLists.txt +++ b/Modules/DiffusionCore/CMakeLists.txt @@ -1,28 +1,28 @@ # With apple gcc 4.2.1 the following waring leads to an build error if boost is enabled if(APPLE) mitkFunctionCheckCAndCXXCompilerFlags("-Wno-error=empty-body" CMAKE_C_FLAGS CMAKE_CXX_FLAGS) endif() MITK_CREATE_MODULE( SUBPROJECTS MITK-Diffusion INCLUDE_DIRS Algorithms Algorithms/Reconstruction Algorithms/Registration Algorithms/Reconstruction/MultishellProcessing Algorithms/Reconstruction/FittingFunctions DicomImport IODataStructures/DiffusionWeightedImages IODataStructures/Properties IODataStructures Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS MitkMapperExt MitkPlanarFigure MitkImageExtraction MitkDICOM MitkMatchPointRegistration PACKAGE_DEPENDS - PUBLIC VTK|FiltersProgrammable Vigra HDF5 + PUBLIC VTK|FiltersProgrammable+IOImage+CommonComputationalGeometry Vigra HDF5 ) if(TARGET ${MODULE_TARGET} AND MITK_USE_OpenMP) target_link_libraries(${MODULE_TARGET} PUBLIC OpenMP::OpenMP_CXX) endif() if(MSVC) mitkFunctionCheckCAndCXXCompilerFlags("/wd4005" CMAKE_C_FLAGS CMAKE_CXX_FLAGS) endif() add_subdirectory(Testing) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/mitkDiffusionImagingConfigure.h.in ${CMAKE_CURRENT_BINARY_DIR}/mitkDiffusionImagingConfigure.h) mitkFunctionGetVersion(${CMAKE_CURRENT_SOURCE_DIR} MITKDIFFUSION) mitkFunctionGetVersionDescription(${CMAKE_CURRENT_SOURCE_DIR} MITKDIFFUSION) configure_file(mitkDiffusionVersion.h.in ${MITK_BINARY_DIR}/mitkDiffusionVersion.h) diff --git a/Modules/DiffusionCore/Rendering/mitkFiberBundleMapper2D.cpp b/Modules/DiffusionCore/Rendering/mitkFiberBundleMapper2D.cpp index bf11bf2..0ed3277 100644 --- a/Modules/DiffusionCore/Rendering/mitkFiberBundleMapper2D.cpp +++ b/Modules/DiffusionCore/Rendering/mitkFiberBundleMapper2D.cpp @@ -1,280 +1,281 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBundleMapper2D.h" #include "mitkBaseRenderer.h" #include "mitkDataNode.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include +#include class vtkShaderCallback : public vtkCommand { public: static vtkShaderCallback *New() { return new vtkShaderCallback; } mitk::BaseRenderer *renderer; mitk::DataNode *node; void Execute(vtkObject *, unsigned long, void*cbo) override { vtkShaderProgram *program = reinterpret_cast(cbo); float fiberOpacity; bool fiberFading = false; float fiberThickness = 0.0; node->GetOpacity(fiberOpacity, nullptr); node->GetFloatProperty("Fiber2DSliceThickness", fiberThickness); node->GetBoolProperty("Fiber2DfadeEFX", fiberFading); program->SetUniformf("fiberOpacity", fiberOpacity); program->SetUniformi("fiberFadingON", fiberFading); program->SetUniformf("fiberThickness", fiberThickness); if (this->renderer) { //get information about current position of views mitk::SliceNavigationController::Pointer sliceContr = renderer->GetSliceNavigationController(); mitk::PlaneGeometry::ConstPointer planeGeo = sliceContr->GetCurrentPlaneGeometry(); //generate according cutting planes based on the view position float planeNormal[3]; planeNormal[0] = planeGeo->GetNormal()[0]; planeNormal[1] = planeGeo->GetNormal()[1]; planeNormal[2] = planeGeo->GetNormal()[2]; float tmp1 = planeGeo->GetOrigin()[0] * planeNormal[0]; float tmp2 = planeGeo->GetOrigin()[1] * planeNormal[1]; float tmp3 = planeGeo->GetOrigin()[2] * planeNormal[2]; float thickness = tmp1 + tmp2 + tmp3; //attention, correct normalvector float a[4]; for (int i = 0; i < 3; ++i) a[i] = planeNormal[i]; a[3] = thickness; program->SetUniform4f("slicingPlane", a); } } vtkShaderCallback() { this->renderer = nullptr; } }; mitk::FiberBundleMapper2D::FiberBundleMapper2D() : m_LineWidth(1) { m_lut = vtkSmartPointer::New(); m_lut->Build(); } mitk::FiberBundleMapper2D::~FiberBundleMapper2D() { } mitk::FiberBundle* mitk::FiberBundleMapper2D::GetInput() { return dynamic_cast< mitk::FiberBundle * > ( GetDataNode()->GetData() ); } void mitk::FiberBundleMapper2D::Update(mitk::BaseRenderer * renderer) { bool visible = true; GetDataNode()->GetVisibility(visible, renderer, "visible"); if ( !visible ) return; // Calculate time step of the input data for the specified renderer (integer value) // this method is implemented in mitkMapper this->CalculateTimeStep( renderer ); //check if updates occured in the node or on the display FBXLocalStorage *localStorage = m_LocalStorageHandler.GetLocalStorage(renderer); //set renderer independent shader properties const DataNode::Pointer node = this->GetDataNode(); float thickness = 2.0; if(!this->GetDataNode()->GetPropertyValue("Fiber2DSliceThickness",thickness)) MITK_INFO << "FIBER2D SLICE THICKNESS PROPERTY ERROR"; bool fiberfading = false; if(!this->GetDataNode()->GetPropertyValue("Fiber2DfadeEFX",fiberfading)) MITK_INFO << "FIBER2D SLICE FADE EFX PROPERTY ERROR"; mitk::FiberBundle* fiberBundle = this->GetInput(); if (fiberBundle==nullptr) return; int lineWidth = 1.0; node->GetIntProperty("LineWidth", lineWidth); if (m_LineWidth!=lineWidth) { m_LineWidth = lineWidth; fiberBundle->RequestUpdate2D(); } vtkProperty *property = localStorage->m_Actor->GetProperty(); property->SetLighting(false); if ( localStorage->m_LastUpdateTimeGetCurrentWorldPlaneGeometryUpdateTime() || localStorage->m_LastUpdateTimeGetUpdateTime2D() ) { this->UpdateShaderParameter(renderer); this->GenerateDataForRenderer( renderer ); } } void mitk::FiberBundleMapper2D::UpdateShaderParameter(mitk::BaseRenderer *) { // see new vtkShaderCallback } // vtkActors and Mappers are feeded here void mitk::FiberBundleMapper2D::GenerateDataForRenderer(mitk::BaseRenderer *renderer) { mitk::FiberBundle* fiberBundle = this->GetInput(); //the handler of local storage gets feeded in this method with requested data for related renderwindow FBXLocalStorage *localStorage = m_LocalStorageHandler.GetLocalStorage(renderer); mitk::DataNode* node = this->GetDataNode(); if (node == nullptr) return; vtkSmartPointer fiberPolyData = fiberBundle->GetFiberPolyData(); if (fiberPolyData == nullptr) return; fiberPolyData->GetPointData()->AddArray(fiberBundle->GetFiberColors()); localStorage->m_Mapper->ScalarVisibilityOn(); localStorage->m_Mapper->SetScalarModeToUsePointFieldData(); localStorage->m_Mapper->SetLookupTable(m_lut); //apply the properties after the slice was set localStorage->m_Actor->GetProperty()->SetOpacity(0.999); localStorage->m_Mapper->SelectColorArray("FIBER_COLORS"); localStorage->m_Mapper->SetInputData(fiberPolyData); - localStorage->m_Mapper->SetVertexShaderCode( + localStorage->m_Actor->GetShaderProperty()->SetVertexShaderCode( "//VTK::System::Dec\n" "attribute vec4 vertexMC;\n" "//VTK::Normal::Dec\n" "uniform mat4 MCDCMatrix;\n" "//VTK::Color::Dec\n" "varying vec4 positionWorld;\n" "varying vec4 colorVertex;\n" "void main(void)\n" "{\n" " colorVertex = scalarColor;\n" " positionWorld = vertexMC;\n" " gl_Position = MCDCMatrix * vertexMC;\n" "}\n" ); - localStorage->m_Mapper->SetFragmentShaderCode( + localStorage->m_Actor->GetShaderProperty()->SetFragmentShaderCode( "//VTK::System::Dec\n" // always start with this line "//VTK::Output::Dec\n" // always have this line in your FS "uniform vec4 slicingPlane;\n" "uniform float fiberThickness;\n" "uniform int fiberFadingON;\n" "uniform float fiberOpacity;\n" "varying vec4 positionWorld;\n" "varying vec4 colorVertex;\n" "out vec4 out_Color;\n" "void main(void)\n" "{\n" " float r1 = dot(positionWorld.xyz, slicingPlane.xyz) - slicingPlane.w;\n" " if (abs(r1) >= fiberThickness)\n" " discard;\n" " if (fiberFadingON != 0)\n" " {\n" " float x = (r1 + fiberThickness) / (fiberThickness*2.0);\n" " x = 1.0 - x;\n" " out_Color = vec4(colorVertex.xyz*x, fiberOpacity);\n" " }\n" " else{\n" " out_Color = vec4(colorVertex.xyz, fiberOpacity);\n" " }\n" "}\n" ); vtkSmartPointer myCallback = vtkSmartPointer::New(); myCallback->renderer = renderer; myCallback->node = this->GetDataNode(); localStorage->m_Mapper->AddObserver(vtkCommand::UpdateShaderEvent,myCallback); localStorage->m_Actor->SetMapper(localStorage->m_Mapper); localStorage->m_Actor->GetProperty()->SetLineWidth(m_LineWidth); // We have been modified => save this for next Update() localStorage->m_LastUpdateTime.Modified(); } vtkProp* mitk::FiberBundleMapper2D::GetVtkProp(mitk::BaseRenderer *renderer) { this->Update(renderer); return m_LocalStorageHandler.GetLocalStorage(renderer)->m_Actor; } void mitk::FiberBundleMapper2D::SetDefaultProperties(mitk::DataNode* node, mitk::BaseRenderer* renderer, bool overwrite) { Superclass::SetDefaultProperties(node, renderer, overwrite); // node->SetProperty("shader",mitk::ShaderProperty::New("mitkShaderFiberClipping")); //add other parameters to propertylist node->AddProperty( "Fiber2DSliceThickness", mitk::FloatProperty::New(1.0f), renderer, overwrite ); node->AddProperty( "Fiber2DfadeEFX", mitk::BoolProperty::New(true), renderer, overwrite ); node->AddProperty( "color", mitk::ColorProperty::New(1.0,1.0,1.0), renderer, overwrite); } mitk::FiberBundleMapper2D::FBXLocalStorage::FBXLocalStorage() { m_Actor = vtkSmartPointer::New(); m_Mapper = vtkSmartPointer::New(); } diff --git a/Modules/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp b/Modules/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp index 935f41a..6592804 100644 --- a/Modules/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp +++ b/Modules/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp @@ -1,183 +1,185 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkIOUtil.h" #include "mitkTestingMacros.h" #include "mitkTestFixture.h" #include #include #include #include #include +#include class mitkDiffusionPropertySerializerTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkDiffusionPropertySerializerTestSuite); MITK_TEST(Equal_SerializeandDeserialize_ReturnsTrue); //MITK_TEST(Equal_DifferentChannels_ReturnFalse); CPPUNIT_TEST_SUITE_END(); private: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ mitk::PropertyList::Pointer propList; //represet image propertylist mitk::BValueMapProperty::Pointer bvaluemap_prop; mitk::GradientDirectionsProperty::Pointer gradientdirection_prop; mitk::MeasurementFrameProperty::Pointer measurementframe_prop; public: /** * @brief Setup Always call this method before each Test-case to ensure correct and new intialization of the used members for a new test case. (If the members are not used in a test, the method does not need to be called). */ void setUp() override { propList = mitk::PropertyList::New(); mitk::BValueMapProperty::BValueMap map; std::vector indices1; indices1.push_back(1); indices1.push_back(2); indices1.push_back(3); indices1.push_back(4); map[0] = indices1; std::vector indices2; indices2.push_back(4); indices2.push_back(3); indices2.push_back(2); indices2.push_back(1); map[1000] = indices2; bvaluemap_prop = mitk::BValueMapProperty::New(map).GetPointer(); propList->SetProperty(mitk::DiffusionPropertyHelper::GetBvaluePropertyName().c_str(), bvaluemap_prop); mitk::GradientDirectionsProperty::GradientDirectionsContainerType::Pointer gdc; gdc = mitk::GradientDirectionsProperty::GradientDirectionsContainerType::New(); double a[3] = {3.0,4.0,1.4}; vnl_vector_fixed vec1; vec1.set(a); gdc->push_back(vec1); double b[3] = {1.0,5.0,123.4}; vnl_vector_fixed vec2; vec2.set(b); gdc->push_back(vec2); double c[3] = {13.0,84.02,13.4}; vnl_vector_fixed vec3; vec3.set(c); gdc->push_back(vec3); gradientdirection_prop = mitk::GradientDirectionsProperty::New(gdc).GetPointer(); propList->ReplaceProperty(mitk::DiffusionPropertyHelper::GetGradientContainerPropertyName().c_str(), gradientdirection_prop); mitk::MeasurementFrameProperty::MeasurementFrameType mft; double row0[3] = {1,0,0}; double row1[3] = {0,1,0}; double row2[3] = {0,0,1}; mft.set_row(0,row0); mft.set_row(1,row1); mft.set_row(2,row2); measurementframe_prop = mitk::MeasurementFrameProperty::New(mft).GetPointer(); propList->ReplaceProperty(mitk::DiffusionPropertyHelper::GetGradientContainerPropertyName().c_str(), measurementframe_prop); } void tearDown() override { } void Equal_SerializeandDeserialize_ReturnsTrue() { assert(propList); /* try to serialize each property in the list, then deserialize again and check for equality */ for (mitk::PropertyList::PropertyMap::const_iterator it = propList->GetMap()->begin(); it != propList->GetMap()->end(); ++it) { const mitk::BaseProperty* prop = it->second; // construct name of serializer class std::string serializername = std::string(prop->GetNameOfClass()) + "Serializer"; std::list allSerializers = itk::ObjectFactoryBase::CreateAllInstance(serializername.c_str()); MITK_TEST_CONDITION(allSerializers.size() > 0, std::string("Creating serializers for ") + serializername); if (allSerializers.size() == 0) { MITK_TEST_OUTPUT( << "serialization not possible, skipping " << prop->GetNameOfClass()); continue; } if (allSerializers.size() > 1) { MITK_TEST_OUTPUT (<< "Warning: " << allSerializers.size() << " serializers found for " << prop->GetNameOfClass() << "testing only the first one."); } mitk::BasePropertySerializer* serializer = dynamic_cast( allSerializers.begin()->GetPointer()); MITK_TEST_CONDITION(serializer != nullptr, serializername + std::string(" is valid")); if (serializer != nullptr) { serializer->SetProperty(prop); - TiXmlElement* valueelement = nullptr; + tinyxml2::XMLElement* valueelement = nullptr; try { - valueelement = serializer->Serialize(); + tinyxml2::XMLDocument doc; + valueelement = serializer->Serialize(doc); // TiXmlPrinter p; // valueelement->Accept(&p); // MITK_INFO << p.CStr(); } catch (...) { } MITK_TEST_CONDITION(valueelement != nullptr, std::string("Serialize property with ") + serializername); if (valueelement == nullptr) { MITK_TEST_OUTPUT( << "serialization failed, skipping deserialization"); continue; } mitk::BaseProperty::Pointer deserializedProp = serializer->Deserialize( valueelement ); MITK_TEST_CONDITION(deserializedProp.IsNotNull(), "serializer created valid property"); if (deserializedProp.IsNotNull()) { MITK_TEST_CONDITION(*(deserializedProp.GetPointer()) == *prop, "deserialized property equals initial property for type " << prop->GetNameOfClass()); } } else { MITK_TEST_OUTPUT( << "created serializer object is of class " << allSerializers.begin()->GetPointer()->GetNameOfClass()) } } // for all properties } }; MITK_TEST_SUITE_REGISTRATION(mitkDiffusionPropertySerializer) diff --git a/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkReader.cpp index 5c3c62f..4ba6eb0 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkReader.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkReader.cpp @@ -1,228 +1,222 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkConnectomicsNetworkReader.h" #include "mitkConnectomicsNetworkDefinitions.h" -#include +#include #include "itksys/SystemTools.hxx" #include #include "mitkGeometry3D.h" #include #include "mitkDiffusionIOMimeTypes.h" +#include namespace mitk { ConnectomicsNetworkReader::ConnectomicsNetworkReader(const ConnectomicsNetworkReader& other) : mitk::AbstractFileReader(other) { } ConnectomicsNetworkReader::ConnectomicsNetworkReader() : mitk::AbstractFileReader( CustomMimeType( mitk::DiffusionIOMimeTypes::CONNECTOMICS_MIMETYPE() ), mitk::DiffusionIOMimeTypes::CONNECTOMICS_MIMETYPE_DESCRIPTION() ) { m_ServiceReg = this->RegisterService(); } ConnectomicsNetworkReader::~ConnectomicsNetworkReader() { } std::vector > ConnectomicsNetworkReader::DoRead() { std::vector > result; std::string location = GetInputLocation(); std::string ext = itksys::SystemTools::GetFilenameLastExtension(location); ext = itksys::SystemTools::LowerCase(ext); if ( location == "") { MITK_ERROR << "No file name specified."; } else if (ext == ".cnf") { try { mitk::ConnectomicsNetwork::Pointer outputNetwork = mitk::ConnectomicsNetwork::New(); - TiXmlDocument doc( location ); - bool loadOkay = doc.LoadFile(); + tinyxml2::XMLDocument doc; + bool loadOkay = doc.LoadFile(location.c_str()); if(!loadOkay) { mitkThrow() << "Could not open file " << location << " for reading."; } - TiXmlHandle hDoc(&doc); - TiXmlElement* pElem; - TiXmlHandle hRoot(nullptr); + tinyxml2::XMLHandle hDoc(&doc); + tinyxml2::XMLElement* pElem; + tinyxml2::XMLHandle hRoot = hDoc.FirstChildElement(); - pElem = hDoc.FirstChildElement().Element(); + pElem = hRoot.ToElement(); // save this for later - hRoot = TiXmlHandle(pElem); + hRoot = tinyxml2::XMLHandle(pElem); //get file version - std::string version(""); - hRoot.Element()->QueryStringAttribute(mitk::ConnectomicsNetworkDefinitions::XML_FILE_VERSION, &version); + std::string version = std::string(pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_FILE_VERSION)); // read geometry - pElem = hRoot.FirstChildElement(mitk::ConnectomicsNetworkDefinitions::XML_GEOMETRY).Element(); + pElem = hRoot.FirstChildElement(mitk::ConnectomicsNetworkDefinitions::XML_GEOMETRY).ToElement(); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); // read origin mitk::Point3D origin; - double temp = 0; - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_X, &temp); - origin[0] = temp; - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Y, &temp); - origin[1] = temp; - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Z, &temp); - origin[2] = temp; + origin[0] = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_X); + origin[1] = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Y); + origin[2] = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Z); geometry->SetOrigin(origin); // read spacing ScalarType spacing[3]; - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_X, &temp); - spacing[0] = temp; - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Y, &temp); - spacing[1] = temp; - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Z, &temp); - spacing[2] = temp; + spacing[0] = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_X); + spacing[1] = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Y); + spacing[2] = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Z); geometry->SetSpacing(spacing); // read transform + double temp = 0; vtkMatrix4x4* m = vtkMatrix4x4::New(); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XX, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XX); m->SetElement(0,0,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XY, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XY); m->SetElement(1,0,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XZ, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XZ); m->SetElement(2,0,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YX, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YX); m->SetElement(0,1,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YY, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YY); m->SetElement(1,1,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YZ, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YZ); m->SetElement(2,1,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZX, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZX); m->SetElement(0,2,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZY, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZY); m->SetElement(1,2,temp); - pElem->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZZ, &temp); + temp = pElem->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZZ); m->SetElement(2,2,temp); m->SetElement(0,3,origin[0]); m->SetElement(1,3,origin[1]); m->SetElement(2,3,origin[2]); m->SetElement(3,3,1); geometry->SetIndexToWorldTransformByVtkMatrix(m); geometry->SetImageGeometry(true); outputNetwork->SetGeometry(geometry); // read network std::map< int, mitk::ConnectomicsNetwork::VertexDescriptorType > idToVertexMap; // read vertices - pElem = hRoot.FirstChildElement(mitk::ConnectomicsNetworkDefinitions::XML_VERTICES).Element(); + pElem = hRoot.FirstChildElement(mitk::ConnectomicsNetworkDefinitions::XML_VERTICES).ToElement(); { // walk through the vertices - TiXmlElement* vertexElement = pElem->FirstChildElement(); + tinyxml2::XMLElement* vertexElement = pElem->FirstChildElement(); for( ; vertexElement; vertexElement=vertexElement->NextSiblingElement()) { std::vector< float > pos; std::string label; int vertexID(0); - vertexElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_X, &temp); + temp = vertexElement->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_X); pos.push_back(temp); - vertexElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Y, &temp); + temp = vertexElement->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Y); pos.push_back(temp); - vertexElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Z, &temp); + temp = vertexElement->DoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Z); pos.push_back(temp); - vertexElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_ID, &vertexID); - vertexElement->QueryStringAttribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_LABEL, &label); + vertexID = vertexElement->IntAttribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_ID); + label = std::string(vertexElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_LABEL)); mitk::ConnectomicsNetwork::VertexDescriptorType newVertex = outputNetwork->AddVertex( vertexID ); outputNetwork->SetLabel( newVertex, label ); outputNetwork->SetCoordinates( newVertex, pos ); if ( idToVertexMap.count( vertexID ) > 0 ) { MITK_ERROR << "Aborting network creation, duplicate vertex ID in file."; return result; } idToVertexMap.insert( std::pair< int, mitk::ConnectomicsNetwork::VertexDescriptorType >( vertexID, newVertex) ); } } // read edges - pElem = hRoot.FirstChildElement(mitk::ConnectomicsNetworkDefinitions::XML_EDGES).Element(); + pElem = hRoot.FirstChildElement(mitk::ConnectomicsNetworkDefinitions::XML_EDGES).ToElement(); { // walk through the edges - TiXmlElement* edgeElement = pElem->FirstChildElement(); + auto edgeElement = pElem->FirstChildElement(); for( ; edgeElement; edgeElement=edgeElement->NextSiblingElement()) { - int edgeID(0), edgeSourceID(0), edgeTargetID(0), edgeWeight(0); + int edgeSourceID(0), edgeTargetID(0), edgeWeight(0); double edgeDoubleWeight(0.0); - edgeElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_ID, &edgeID); - edgeElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_SOURCE_ID, &edgeSourceID); - edgeElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_TARGET_ID, &edgeTargetID); - edgeElement->Attribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_FIBERCOUNT_ID, &edgeWeight); + edgeSourceID = edgeElement->IntAttribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_SOURCE_ID); + edgeTargetID = edgeElement->IntAttribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_TARGET_ID); + edgeWeight = edgeElement->IntAttribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_FIBERCOUNT_ID); + if(version == "0.1") { // in version 0.1 only the int weight was saved, assume double weight to be one by default edgeDoubleWeight = 1.0; } else { edgeElement->QueryDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_EDGE_DOUBLE_WEIGHT_ID, &edgeDoubleWeight); } mitk::ConnectomicsNetwork::VertexDescriptorType source = idToVertexMap.find( edgeSourceID )->second; mitk::ConnectomicsNetwork::VertexDescriptorType target = idToVertexMap.find( edgeTargetID )->second; outputNetwork->AddEdge( source, target, edgeSourceID, edgeTargetID, edgeWeight, edgeDoubleWeight); } } outputNetwork->UpdateBounds(); result.push_back(outputNetwork.GetPointer()); MITK_INFO << "Network read"; } catch (mitk::Exception& e) { MITK_ERROR << e.GetDescription(); } catch(...) { MITK_ERROR << "Unknown error occured while trying to read file."; } } return result; } } //namespace MITK mitk::ConnectomicsNetworkReader* mitk::ConnectomicsNetworkReader::Clone() const { return new ConnectomicsNetworkReader(*this); } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkWriter.cpp b/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkWriter.cpp index abc33fa..a7a07b9 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkWriter.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkConnectomicsNetworkWriter.cpp @@ -1,144 +1,145 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkConnectomicsNetworkWriter.h" #include "mitkConnectomicsNetworkDefinitions.h" -#include +#include #include "itksys/SystemTools.hxx" #include "mitkDiffusionIOMimeTypes.h" mitk::ConnectomicsNetworkWriter::ConnectomicsNetworkWriter() : AbstractFileWriter(mitk::ConnectomicsNetwork::GetStaticNameOfClass(), CustomMimeType( mitk::DiffusionIOMimeTypes::CONNECTOMICS_MIMETYPE() ), mitk::DiffusionIOMimeTypes::CONNECTOMICS_MIMETYPE_DESCRIPTION() ) { RegisterService(); } mitk::ConnectomicsNetworkWriter::ConnectomicsNetworkWriter(const mitk::ConnectomicsNetworkWriter& other) : AbstractFileWriter(other) { } mitk::ConnectomicsNetworkWriter::~ConnectomicsNetworkWriter() {} mitk::ConnectomicsNetworkWriter* mitk::ConnectomicsNetworkWriter::Clone() const { return new ConnectomicsNetworkWriter(*this); } void mitk::ConnectomicsNetworkWriter::Write() { MITK_INFO << "Writing connectomics network"; InputType::ConstPointer input = dynamic_cast(this->GetInput()); if (input.IsNull() ) { MITK_ERROR <<"Sorry, input to ConnectomicsNetworkWriter is nullptr!"; return; } if ( this->GetOutputLocation().empty() ) { MITK_ERROR << "Sorry, filename has not been set!" ; return ; } std::string ext = itksys::SystemTools::GetFilenameLastExtension(this->GetOutputLocation()); ext = itksys::SystemTools::LowerCase(ext); // default extension is .cnf if(ext == "") { ext = ".cnf"; this->SetOutputLocation(this->GetOutputLocation() + ext); } if (ext == ".cnf") { // Get geometry of the network mitk::BaseGeometry* geometry = input->GetGeometry(); // Create XML document - TiXmlDocument documentXML; + tinyxml2::XMLDocument documentXML; { // begin document - TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); // TODO what to write here? encoding? etc.... - documentXML.LinkEndChild( declXML ); + //tinyxml2::XMLDeclaration* declXML = new tinyxml2::XMLDeclaration( "1.0", "", "" ); // TODO what to write here? encoding? etc.... + //documentXML.LinkEndChild( declXML ); - TiXmlElement* mainXML = new TiXmlElement(mitk::ConnectomicsNetworkDefinitions::XML_CONNECTOMICS_FILE); + auto mainXML = documentXML.NewElement(mitk::ConnectomicsNetworkDefinitions::XML_CONNECTOMICS_FILE); mainXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_FILE_VERSION, mitk::ConnectomicsNetworkDefinitions::VERSION_STRING); - documentXML.LinkEndChild(mainXML); + //documentXML.LinkEndChild(mainXML); + documentXML.InsertFirstChild(mainXML); - TiXmlElement* geometryXML = new TiXmlElement(mitk::ConnectomicsNetworkDefinitions::XML_GEOMETRY); + auto geometryXML = documentXML.NewElement(mitk::ConnectomicsNetworkDefinitions::XML_GEOMETRY); { // begin geometry - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XX, geometry->GetMatrixColumn(0)[0]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XY, geometry->GetMatrixColumn(0)[1]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XZ, geometry->GetMatrixColumn(0)[2]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YX, geometry->GetMatrixColumn(1)[0]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YY, geometry->GetMatrixColumn(1)[1]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YZ, geometry->GetMatrixColumn(1)[2]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZX, geometry->GetMatrixColumn(2)[0]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZY, geometry->GetMatrixColumn(2)[1]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZZ, geometry->GetMatrixColumn(2)[2]); - - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_X, geometry->GetOrigin()[0]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Y, geometry->GetOrigin()[1]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Z, geometry->GetOrigin()[2]); - - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_X, geometry->GetSpacing()[0]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Y, geometry->GetSpacing()[1]); - geometryXML->SetDoubleAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Z, geometry->GetSpacing()[2]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XX, geometry->GetMatrixColumn(0)[0]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XY, geometry->GetMatrixColumn(0)[1]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_XZ, geometry->GetMatrixColumn(0)[2]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YX, geometry->GetMatrixColumn(1)[0]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YY, geometry->GetMatrixColumn(1)[1]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_YZ, geometry->GetMatrixColumn(1)[2]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZX, geometry->GetMatrixColumn(2)[0]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZY, geometry->GetMatrixColumn(2)[1]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_MATRIX_ZZ, geometry->GetMatrixColumn(2)[2]); + + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_X, geometry->GetOrigin()[0]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Y, geometry->GetOrigin()[1]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_ORIGIN_Z, geometry->GetOrigin()[2]); + + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_X, geometry->GetSpacing()[0]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Y, geometry->GetSpacing()[1]); + geometryXML->SetAttribute(mitk::ConnectomicsNetworkDefinitions::XML_SPACING_Z, geometry->GetSpacing()[2]); } // end geometry - mainXML->LinkEndChild(geometryXML); + mainXML->InsertEndChild(geometryXML); - TiXmlElement* verticesXML = new TiXmlElement(mitk::ConnectomicsNetworkDefinitions::XML_VERTICES); + auto verticesXML = documentXML.NewElement(mitk::ConnectomicsNetworkDefinitions::XML_VERTICES); { // begin vertices section VertexVectorType vertexVector = dynamic_cast(this->GetInput())->GetVectorOfAllNodes(); for( unsigned int index = 0; index < vertexVector.size(); index++ ) { // not localized as of yet TODO - TiXmlElement* vertexXML = new TiXmlElement(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX ); + auto vertexXML = documentXML.NewElement(mitk::ConnectomicsNetworkDefinitions::XML_VERTEX ); vertexXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_ID , vertexVector[ index ].id ); - vertexXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_LABEL , vertexVector[ index ].label ); - vertexXML->SetDoubleAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_X , vertexVector[ index ].coordinates[0] ); - vertexXML->SetDoubleAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Y , vertexVector[ index ].coordinates[1] ); - vertexXML->SetDoubleAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Z , vertexVector[ index ].coordinates[2] ); - verticesXML->LinkEndChild(vertexXML); + vertexXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_LABEL , vertexVector[ index ].label.c_str() ); + vertexXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_X , vertexVector[ index ].coordinates[0] ); + vertexXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Y , vertexVector[ index ].coordinates[1] ); + vertexXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_VERTEX_Z , vertexVector[ index ].coordinates[2] ); + verticesXML->InsertEndChild(vertexXML); } } // end vertices section - mainXML->LinkEndChild(verticesXML); + mainXML->InsertEndChild(verticesXML); - TiXmlElement* edgesXML = new TiXmlElement(mitk::ConnectomicsNetworkDefinitions::XML_EDGES); + auto edgesXML = documentXML.NewElement(mitk::ConnectomicsNetworkDefinitions::XML_EDGES); { // begin edges section EdgeVectorType edgeVector = dynamic_cast(this->GetInput())->GetVectorOfAllEdges(); for(unsigned int index = 0; index < edgeVector.size(); index++ ) { - TiXmlElement* edgeXML = new TiXmlElement(mitk::ConnectomicsNetworkDefinitions::XML_EDGE ); + auto edgeXML = documentXML.NewElement(mitk::ConnectomicsNetworkDefinitions::XML_EDGE ); edgeXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_EDGE_ID , index ); edgeXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_EDGE_SOURCE_ID , edgeVector[ index ].second.sourceId ); edgeXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_EDGE_TARGET_ID , edgeVector[ index ].second.targetId ); edgeXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_EDGE_FIBERCOUNT_ID , edgeVector[ index ].second.fiber_count ); - edgeXML->SetDoubleAttribute( mitk::ConnectomicsNetworkDefinitions::XML_EDGE_DOUBLE_WEIGHT_ID , edgeVector[ index ].second.edge_weight ); - edgesXML->LinkEndChild(edgeXML); + edgeXML->SetAttribute( mitk::ConnectomicsNetworkDefinitions::XML_EDGE_DOUBLE_WEIGHT_ID , edgeVector[ index ].second.edge_weight ); + edgesXML->InsertEndChild(edgeXML); } } // end edges section - mainXML->LinkEndChild(edgesXML); + mainXML->InsertEndChild(edgesXML); } // end document documentXML.SaveFile( this->GetOutputLocation().c_str() ); MITK_INFO << "Connectomics network written"; } } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleDicomReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleDicomReader.cpp index 213a204..6987b09 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleDicomReader.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleDicomReader.cpp @@ -1,166 +1,165 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBundleDicomReader.h" #include #include #include #include #include #include #include #include #include #include #include #include -#include #include #include #include #include "mitkDiffusionIOMimeTypes.h" #include #include #include "dcmtk/config/osconfig.h" /* make sure OS specific configuration is included first */ #include "dcmtk/dcmtract/trctractographyresults.h" #include "dcmtk/dcmtract/trctrack.h" mitk::FiberBundleDicomReader::FiberBundleDicomReader() : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_DICOM_MIMETYPE_NAME(), "DICOM Fiber Bundle Reader" ) { m_ServiceReg = this->RegisterService(); } mitk::FiberBundleDicomReader::FiberBundleDicomReader(const FiberBundleDicomReader &other) :mitk::AbstractFileReader(other) { } mitk::FiberBundleDicomReader * mitk::FiberBundleDicomReader::Clone() const { return new FiberBundleDicomReader(*this); } std::vector > mitk::FiberBundleDicomReader::DoRead() { std::vector > output_fibs; try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::string filename = this->GetInputLocation(); OFCondition result; TrcTractographyResults *trc = nullptr; result = TrcTractographyResults::loadFile(filename.c_str(), trc); if (result.bad()) mitkThrow() << "Unable to load tractography dicom file: " << result.text(); OFString val = "-"; trc->getPatient().getPatientName(val); MITK_INFO << "Patient Name: " << val; val = "-"; trc->getStudy().getStudyInstanceUID(val); MITK_INFO << "Study : " << val; val = "-"; trc->getSeries().getSeriesInstanceUID(val); MITK_INFO << "Series : " << val; val = "-"; trc->getSOPCommon().getSOPInstanceUID(val); MITK_INFO << "Instance : " << val; val = "-"; MITK_INFO << "-------------------------------------------------------------------------"; size_t numTrackSets = trc->getNumberOfTrackSets(); OFVector& sets = trc->getTrackSets(); for (size_t ts = 0; ts < numTrackSets; ts++) { size_t numTracks = sets[ts]->getNumberOfTracks(); MITK_INFO << " Track Set #" << ts << ": " << numTracks << " Tracks, " << sets[ts]->getNumberOfTrackSetStatistics() << " Track Set Statistics, " << sets[ts]->getNumberOfTrackStatistics() << " Track Statistics, " << sets[ts]->getNumberOfMeasurements() << " Measurements "; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (size_t t = 0; t < numTracks; t++) { vtkSmartPointer container = vtkSmartPointer::New(); TrcTrack* track = sets[ts]->getTracks()[t]; const Float32* vals = nullptr; size_t numPoints = track->getTrackData(vals); for (size_t v = 0; v < numPoints; ++v) { vtkIdType id = vtkNewPoints->InsertNextPoint(vals[v*3],vals[v*3+1],vals[v*3+2]); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(vtkNewPoints); fiberPolyData->SetLines(vtkNewCells); FiberBundle::Pointer fib = FiberBundle::New(fiberPolyData); CodeSequenceMacro algoCode = sets[ts]->getTrackingAlgorithmIdentification().at(0)->getAlgorithmFamilyCode(); val = "0"; algoCode.getCodeValue(val); fib->GetPropertyList()->SetStringProperty("DICOM.algo_family_code.value",val.c_str()); val = "0"; algoCode.getCodeMeaning(val); fib->GetPropertyList()->SetStringProperty("DICOM.algo_family_code.meaning",val.c_str()); CodeSequenceMacro modelCode = sets[ts]->getDiffusionModelCode(); val = "0"; modelCode.getCodeValue(val); fib->GetPropertyList()->SetStringProperty("DICOM.model_code.value",val.c_str()); val = "0"; modelCode.getCodeMeaning(val); fib->GetPropertyList()->SetStringProperty("DICOM.model_code.meaning",val.c_str()); CodeWithModifiers anatomy = sets[ts]->getTrackSetAnatomy(); val = "0"; anatomy.getCodeValue(val); fib->GetPropertyList()->SetStringProperty("DICOM.anatomy.value",val.c_str()); val = "0"; anatomy.getCodeMeaning(val); fib->GetPropertyList()->SetStringProperty("DICOM.anatomy.meaning",val.c_str()); val = "0"; trc->getPatient().getPatientID(val); fib->GetPropertyList()->SetStringProperty("DICOM.patient_id",val.c_str()); val = "0"; trc->getPatient().getPatientName(val); fib->GetPropertyList()->SetStringProperty("DICOM.patient_name",val.c_str()); val = "0"; trc->getStudy().getStudyInstanceUID(val); fib->GetPropertyList()->SetStringProperty("DICOM.study_instance_uid",val.c_str()); val = "0"; trc->getSeries().getSeriesInstanceUID(val); fib->GetPropertyList()->SetStringProperty("DICOM.series_instance_uid",val.c_str()); val = "0"; trc->getSOPCommon().getSOPInstanceUID(val); fib->GetPropertyList()->SetStringProperty("DICOM.sop_instance_uid",val.c_str()); val = "0"; trc->getSOPCommon().getSOPClassUID(val); fib->GetPropertyList()->SetStringProperty("DICOM.sop_class_uid",val.c_str()); val = "0"; trc->getFrameOfReference().getFrameOfReferenceUID(val); fib->GetPropertyList()->SetStringProperty("DICOM.frame_of_reference_uid",val.c_str()); output_fibs.push_back(fib.GetPointer()); MITK_INFO << "Fiber bundle read"; } delete trc; setlocale(LC_ALL, currLocale.c_str()); return output_fibs; } catch(...) { throw; } return output_fibs; } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTckReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTckReader.cpp index 6099729..d370b36 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTckReader.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTckReader.cpp @@ -1,164 +1,163 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBundleTckReader.h" #include #include #include #include #include #include #include #include #include #include #include #include -#include #include #include #include #include "mitkDiffusionIOMimeTypes.h" #include #include mitk::FiberBundleTckReader::FiberBundleTckReader() : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_TCK_MIMETYPE_NAME(), "tck Fiber Bundle Reader (MRtrix format)" ) { m_ServiceReg = this->RegisterService(); } mitk::FiberBundleTckReader::FiberBundleTckReader(const FiberBundleTckReader &other) :mitk::AbstractFileReader(other) { } mitk::FiberBundleTckReader * mitk::FiberBundleTckReader::Clone() const { return new FiberBundleTckReader(*this); } std::vector > mitk::FiberBundleTckReader::DoRead() { std::vector > result; try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::string filename = this->GetInputLocation(); std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); if (ext==".tck") { MITK_INFO << "Loading tractogram (MRtrix format): " << itksys::SystemTools::GetFilenameName(filename); std::FILE* filePointer = std::fopen(filename.c_str(),"r+b"); std::string header = ""; bool header_end = false; int count = 0; while (!header_end) { char c; count += std::fread(&c, 1, 1, filePointer); header += c; if (header.size() >= 3 && header.compare(header.size() - 3, 3, "END") == 0) header_end = true; } MITK_INFO << "TCK Header:"; MITK_INFO << header; int header_size = -1; try { std::string delimiter = "file: . "; size_t pos = 0; pos = header.find(delimiter); header.erase(0, pos + delimiter.length()); pos = header.find("\n"); header_size = boost::lexical_cast(header.substr(0, pos)); } catch(...) { } if (header_size==-1) mitkThrow() << "Could not parse header size from " << filename; std::fseek ( filePointer , header_size , SEEK_SET ); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer container = vtkSmartPointer::New(); float tmp[3]; count = 1; int fiber_length = 0; while (std::fread((char*)tmp, 1, 12, filePointer)==12) { if (std::isinf(tmp[0]) || std::isinf(tmp[1]) || std::isinf(tmp[1])) break; else if (std::isnan(tmp[0]) || std::isnan(tmp[1]) || std::isnan(tmp[2])) { count++; fiber_length = 0; vtkNewCells->InsertNextCell(container); container = vtkSmartPointer::New(); } else { fiber_length++; vtkIdType id = vtkNewPoints->InsertNextPoint(tmp); container->GetPointIds()->InsertNextId(id); } } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(vtkNewPoints); fiberPolyData->SetLines(vtkNewCells); // transform polydata from RAS (MRtrix) to LPS (MITK) mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); vtkSmartPointer< vtkMatrix4x4 > matrix = vtkSmartPointer< vtkMatrix4x4 >::New(); matrix->Identity(); matrix->SetElement(0,0,-matrix->GetElement(0,0)); matrix->SetElement(1,1,-matrix->GetElement(1,1)); geometry->SetIndexToWorldTransformByVtkMatrix(matrix); vtkSmartPointer transformFilter = vtkSmartPointer::New(); transformFilter->SetInputData(fiberPolyData); transformFilter->SetTransform(geometry->GetVtkTransform()); transformFilter->Update(); FiberBundle::Pointer fib = FiberBundle::New(transformFilter->GetOutput()); result.push_back(fib.GetPointer()); } setlocale(LC_ALL, currLocale.c_str()); MITK_INFO << "Fiber bundle read"; return result; } catch(...) { throw; } return result; } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp index 671d58b..ea66f0d 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleTrackVisReader.cpp @@ -1,95 +1,94 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBundleTrackVisReader.h" #include #include #include #include #include #include #include #include #include #include #include #include -#include #include #include #include #include "mitkDiffusionIOMimeTypes.h" mitk::FiberBundleTrackVisReader::FiberBundleTrackVisReader() : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_TRK_MIMETYPE_NAME(), "TrackVis Fiber Bundle Reader" ) { Options defaultOptions; defaultOptions["Apply index to world transform stored in the TRK header"] = true; defaultOptions["Print header"] = false; this->SetDefaultOptions(defaultOptions); m_ServiceReg = this->RegisterService(); } mitk::FiberBundleTrackVisReader::FiberBundleTrackVisReader(const FiberBundleTrackVisReader &other) :mitk::AbstractFileReader(other) { } mitk::FiberBundleTrackVisReader * mitk::FiberBundleTrackVisReader::Clone() const { return new FiberBundleTrackVisReader(*this); } std::vector > mitk::FiberBundleTrackVisReader::DoRead() { std::vector > result; try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::string filename = this->GetInputLocation(); MITK_INFO << "Loading tractogram (TrackVis format): " << itksys::SystemTools::GetFilenameName(filename); std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); if (ext==".trk") { Options options = this->GetOptions(); bool apply_matrix = us::any_cast(options["Apply index to world transform stored in the TRK header"]); bool print_header = us::any_cast(options["Print header"]); FiberBundle::Pointer mitk_fib = FiberBundle::New(); TrackVisFiberReader reader; reader.open(this->GetInputLocation().c_str()); reader.read(mitk_fib.GetPointer(), apply_matrix, print_header); result.push_back(mitk_fib.GetPointer()); return result; } setlocale(LC_ALL, currLocale.c_str()); MITK_INFO << "Fiber bundle read"; } catch(...) { throw; } return result; } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkReader.cpp index 15d33bf..1e320a2 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkReader.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkReader.cpp @@ -1,284 +1,170 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBundleVtkReader.h" #include #include #include #include #include #include #include #include #include #include #include #include -#include +#include #include #include #include #include #include #include "mitkDiffusionIOMimeTypes.h" +#include mitk::FiberBundleVtkReader::FiberBundleVtkReader() : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::FIBERBUNDLE_VTK_MIMETYPE_NAME(), "VTK Fiber Bundle Reader" ) { m_ServiceReg = this->RegisterService(); } mitk::FiberBundleVtkReader::FiberBundleVtkReader(const FiberBundleVtkReader &other) :mitk::AbstractFileReader(other) { } mitk::FiberBundleVtkReader * mitk::FiberBundleVtkReader::Clone() const { return new FiberBundleVtkReader(*this); } std::vector > mitk::FiberBundleVtkReader::DoRead() { std::vector > result; const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::string filename = this->GetInputLocation(); std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); try { MITK_INFO << "Loading tractogram (VTK format): " << itksys::SystemTools::GetFilenameName(filename); vtkSmartPointer reader = vtkSmartPointer::New(); reader->SetFileName( this->GetInputLocation().c_str() ); if (reader->IsFilePolyData()) { reader->Update(); if ( reader->GetOutput() != nullptr ) { vtkSmartPointer fiberPolyData = reader->GetOutput(); FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData); vtkSmartPointer weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS")); if (weights!=nullptr) { // float weight=0; for (int i=0; iGetNumberOfValues(); i++) { if (weights->GetValue(i)<0.0) { MITK_ERROR << "Fiber weight<0 detected! Setting value to 0."; weights->SetValue(i,0); } } // if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) // { // MITK_INFO << "Weight: " << weights->GetValue(i); // weight = weights->GetValue(i); // } fiberBundle->SetFiberWeights(weights); } vtkSmartPointer fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS")); if (fiberColors!=nullptr) fiberBundle->SetFiberColors(fiberColors); result.push_back(fiberBundle.GetPointer()); return result; } } else MITK_INFO << "File is not VTK format."; } catch(...) { throw; } try { MITK_INFO << "Trying to load fiber file as VTP format."; vtkSmartPointer reader = vtkSmartPointer::New(); reader->SetFileName( this->GetInputLocation().c_str() ); if ( reader->CanReadFile(this->GetInputLocation().c_str()) ) { reader->Update(); if ( reader->GetOutput() != nullptr ) { vtkSmartPointer fiberPolyData = reader->GetOutput(); FiberBundle::Pointer fiberBundle = FiberBundle::New(fiberPolyData); vtkSmartPointer weights = vtkFloatArray::SafeDownCast(fiberPolyData->GetCellData()->GetArray("FIBER_WEIGHTS")); if (weights!=nullptr) { // float weight=0; // for (int i=0; iGetSize(); i++) // if (!mitk::Equal(weights->GetValue(i),weight,0.00001)) // { // MITK_INFO << "Weight: " << weights->GetValue(i); // weight = weights->GetValue(i); // } fiberBundle->SetFiberWeights(weights); } vtkSmartPointer fiberColors = vtkUnsignedCharArray::SafeDownCast(fiberPolyData->GetPointData()->GetArray("FIBER_COLORS")); if (fiberColors!=nullptr) fiberBundle->SetFiberColors(fiberColors); result.push_back(fiberBundle.GetPointer()); return result; } } else MITK_INFO << "File is not VTP format."; } catch(...) { throw; } - try - { - MITK_INFO << "Trying to load fiber file as deprecated XML format."; - vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); - vtkSmartPointer cellArray = vtkSmartPointer::New(); - vtkSmartPointer points = vtkSmartPointer::New(); - TiXmlDocument doc( this->GetInputLocation().c_str() ); - if(doc.LoadFile()) - { - TiXmlHandle hDoc(&doc); - TiXmlElement* pElem; - TiXmlHandle hRoot(nullptr); - pElem = hDoc.FirstChildElement().Element(); - // save this for later - hRoot = TiXmlHandle(pElem); - pElem = hRoot.FirstChildElement("geometry").Element(); - // read geometry - mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); - // read origin - mitk::Point3D origin; - double temp = 0; - pElem->Attribute("origin_x", &temp); - origin[0] = temp; - pElem->Attribute("origin_y", &temp); - origin[1] = temp; - pElem->Attribute("origin_z", &temp); - origin[2] = temp; - geometry->SetOrigin(origin); - // read spacing - ScalarType spacing[3]; - pElem->Attribute("spacing_x", &temp); - spacing[0] = temp; - pElem->Attribute("spacing_y", &temp); - spacing[1] = temp; - pElem->Attribute("spacing_z", &temp); - spacing[2] = temp; - geometry->SetSpacing(spacing); - // read transform - vtkMatrix4x4* m = vtkMatrix4x4::New(); - pElem->Attribute("xx", &temp); - m->SetElement(0,0,temp); - pElem->Attribute("xy", &temp); - m->SetElement(1,0,temp); - pElem->Attribute("xz", &temp); - m->SetElement(2,0,temp); - pElem->Attribute("yx", &temp); - m->SetElement(0,1,temp); - pElem->Attribute("yy", &temp); - m->SetElement(1,1,temp); - pElem->Attribute("yz", &temp); - m->SetElement(2,1,temp); - pElem->Attribute("zx", &temp); - m->SetElement(0,2,temp); - pElem->Attribute("zy", &temp); - m->SetElement(1,2,temp); - pElem->Attribute("zz", &temp); - m->SetElement(2,2,temp); - m->SetElement(0,3,origin[0]); - m->SetElement(1,3,origin[1]); - m->SetElement(2,3,origin[2]); - m->SetElement(3,3,1); - geometry->SetIndexToWorldTransformByVtkMatrix(m); - // read bounds - float bounds[] = {0, 0, 0, 0, 0, 0}; - pElem->Attribute("size_x", &temp); - bounds[1] = temp; - pElem->Attribute("size_y", &temp); - bounds[3] = temp; - pElem->Attribute("size_z", &temp); - bounds[5] = temp; - geometry->SetFloatBounds(bounds); - geometry->SetImageGeometry(true); - pElem = hRoot.FirstChildElement("fiber_bundle").FirstChild().Element(); - for( ; pElem ; pElem=pElem->NextSiblingElement()) - { - TiXmlElement* pElem2 = pElem->FirstChildElement(); - vtkSmartPointer container = vtkSmartPointer::New(); - for( ; pElem2; pElem2=pElem2->NextSiblingElement()) - { - Point3D point; - pElem2->Attribute("pos_x", &temp); - point[0] = temp; - pElem2->Attribute("pos_y", &temp); - point[1] = temp; - pElem2->Attribute("pos_z", &temp); - point[2] = temp; - geometry->IndexToWorld(point, point); - vtkIdType id = points->InsertNextPoint(point.GetDataPointer()); - container->GetPointIds()->InsertNextId(id); - } - cellArray->InsertNextCell(container); - } - fiberPolyData->SetPoints(points); - fiberPolyData->SetLines(cellArray); - vtkSmartPointer cleaner = vtkSmartPointer::New(); - cleaner->SetInputData(fiberPolyData); - cleaner->Update(); - fiberPolyData = cleaner->GetOutput(); - FiberBundle::Pointer image = FiberBundle::New(fiberPolyData); - result.push_back(image.GetPointer()); - return result; - } - else - { - MITK_INFO << "File is not deprectaed XML format."; - } - - setlocale(LC_ALL, currLocale.c_str()); - MITK_INFO << "Fiber bundle read"; - } - catch(...) - { - throw; - } - throw "Selected file is no vtk readable fiber format (binary or ascii vtk or vtp file)."; return result; } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkWriter.cpp b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkWriter.cpp index c7654c4..ba48a9a 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkWriter.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkFiberBundleVtkWriter.cpp @@ -1,171 +1,172 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBundleVtkWriter.h" #include #include #include #include #include #include #include #include #include #include #include #include "mitkDiffusionIOMimeTypes.h" +#include mitk::FiberBundleVtkWriter::FiberBundleVtkWriter() : mitk::AbstractFileWriter(mitk::FiberBundle::GetStaticNameOfClass(), mitk::DiffusionIOMimeTypes::FIBERBUNDLE_VTK_MIMETYPE_NAME(), "VTK Fiber Bundle Writer") { Options defaultOptions; defaultOptions["Save as binary file"] = true; defaultOptions["Save as xml file (vtp style)"] = false; defaultOptions["Save color information"] = true; defaultOptions["Save fiber weights"] = true; this->SetDefaultOptions(defaultOptions); RegisterService(); } mitk::FiberBundleVtkWriter::FiberBundleVtkWriter(const mitk::FiberBundleVtkWriter & other) :mitk::AbstractFileWriter(other) {} mitk::FiberBundleVtkWriter::~FiberBundleVtkWriter() {} mitk::FiberBundleVtkWriter * mitk::FiberBundleVtkWriter::Clone() const { return new mitk::FiberBundleVtkWriter(*this); } void mitk::FiberBundleVtkWriter::Write() { std::ostream* out; std::ofstream outStream; if( this->GetOutputStream() ) { out = this->GetOutputStream(); }else{ outStream.open( this->GetOutputLocation().c_str() ); out = &outStream; } if ( !out->good() ) { mitkThrow() << "Stream not good."; } try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::locale I("C"); out->imbue(I); std::string filename = this->GetOutputLocation().c_str(); mitk::FiberBundle::ConstPointer input = dynamic_cast(this->GetInput()); std::string ext = itksys::SystemTools::GetFilenameLastExtension(this->GetOutputLocation().c_str()); Options options = this->GetOptions(); vtkSmartPointer weights = nullptr; vtkSmartPointer fiberColors = nullptr; vtkSmartPointer fibPoly = input->GetFiberPolyData(); if (us::any_cast(options["Save fiber weights"])) { MITK_INFO << "Adding fiber weight information"; fibPoly->GetCellData()->AddArray(input->GetFiberWeights()); } else if (fibPoly->GetCellData()->HasArray("FIBER_WEIGHTS")) { weights = input->GetFiberWeights(); fibPoly->GetCellData()->RemoveArray("FIBER_WEIGHTS"); } if (us::any_cast(options["Save color information"])) { MITK_INFO << "Adding color information"; fibPoly->GetPointData()->AddArray(input->GetFiberColors()); } else if (fibPoly->GetPointData()->HasArray("FIBER_COLORS")) { fiberColors = input->GetFiberColors(); fibPoly->GetPointData()->RemoveArray("FIBER_COLORS"); } // default extension is .fib if(ext == "") { ext = ".fib"; this->SetOutputLocation(this->GetOutputLocation() + ext); } if (us::any_cast(options["Save as xml file (vtp style)"])) { vtkSmartPointer writer = vtkSmartPointer::New(); writer->SetInputData(fibPoly); writer->SetFileName(filename.c_str()); if (us::any_cast(options["Save as binary file"])) { MITK_INFO << "Writing fiber bundle as vtk binary file"; writer->SetDataModeToBinary(); } else { MITK_INFO << "Writing fiber bundle as vtk ascii file"; writer->SetDataModeToAscii(); } writer->Write(); } else { vtkSmartPointer writer = vtkSmartPointer::New(); writer->SetInputData(fibPoly); writer->SetFileName(filename.c_str()); if (us::any_cast(options["Save as binary file"])) { MITK_INFO << "Writing fiber bundle as vtk binary file"; writer->SetFileTypeToBinary(); } else { MITK_INFO << "Writing fiber bundle as vtk ascii file"; writer->SetFileTypeToASCII(); } writer->Write(); } // restore arrays if ( weights!=nullptr ) fibPoly->GetPointData()->AddArray(weights); if (fiberColors != nullptr) fibPoly->GetPointData()->AddArray(fiberColors); setlocale(LC_ALL, currLocale.c_str()); MITK_INFO << "VTK Fiber bundle written to " << filename; } catch(...) { throw; } } diff --git a/Modules/DiffusionIO/ReaderWriter/mitkPlanarFigureCompositeReader.cpp b/Modules/DiffusionIO/ReaderWriter/mitkPlanarFigureCompositeReader.cpp index 1240062..365bf93 100644 --- a/Modules/DiffusionIO/ReaderWriter/mitkPlanarFigureCompositeReader.cpp +++ b/Modules/DiffusionIO/ReaderWriter/mitkPlanarFigureCompositeReader.cpp @@ -1,100 +1,99 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPlanarFigureCompositeReader.h" #include #include -#include #include #include #include #include "mitkDiffusionIOMimeTypes.h" #include #define RAPIDXML_NO_EXCEPTIONS #include #include #include #include mitk::PlanarFigureCompositeReader::PlanarFigureCompositeReader() : mitk::AbstractFileReader( mitk::DiffusionIOMimeTypes::PLANARFIGURECOMPOSITE_MIMETYPE(), "Planar Figure Composite Reader" ) { m_ServiceReg = this->RegisterService(); } mitk::PlanarFigureCompositeReader::PlanarFigureCompositeReader(const PlanarFigureCompositeReader &other) :mitk::AbstractFileReader(other) { } mitk::PlanarFigureCompositeReader * mitk::PlanarFigureCompositeReader::Clone() const { return new PlanarFigureCompositeReader(*this); } std::vector > mitk::PlanarFigureCompositeReader::DoRead() { std::vector > result; try { const std::string& locale = "C"; const std::string& currLocale = setlocale( LC_ALL, nullptr ); setlocale(LC_ALL, locale.c_str()); std::string filename = this->GetInputLocation(); std::string ext = itksys::SystemTools::GetFilenameLastExtension(filename); ext = itksys::SystemTools::LowerCase(ext); boost::property_tree::ptree tree; boost::property_tree::xml_parser::read_xml(filename, tree); int comptype = tree.get("comptype"); mitk::PlanarFigureComposite::Pointer pfc = mitk::PlanarFigureComposite::New(); switch(comptype) { case 0: pfc->setOperationType(mitk::PlanarFigureComposite::AND); MITK_INFO << "loading AND composition"; break; case 1: pfc->setOperationType(mitk::PlanarFigureComposite::OR); MITK_INFO << "loading OR composition"; break; case 2: pfc->setOperationType(mitk::PlanarFigureComposite::NOT); MITK_INFO << "loading NOT composition"; break; default: MITK_ERROR << filename << " contains no valid composition type!"; } std::vector > result; result.push_back(pfc.GetPointer()); setlocale(LC_ALL, currLocale.c_str()); return result; } catch(...) { throw; } return result; } diff --git a/Modules/FiberTracking/Algorithms/GibbsTracking/mitkSphereInterpolator.cpp b/Modules/FiberTracking/Algorithms/GibbsTracking/mitkSphereInterpolator.cpp index bec1b21..ce206c4 100644 --- a/Modules/FiberTracking/Algorithms/GibbsTracking/mitkSphereInterpolator.cpp +++ b/Modules/FiberTracking/Algorithms/GibbsTracking/mitkSphereInterpolator.cpp @@ -1,173 +1,173 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkSphereInterpolator.h" #include #include #include #include #include #include #include #include static const std::string BaryCoordsFileName = "FiberTrackingLUTBaryCoords.bin"; static const std::string IndicesFileName = "FiberTrackingLUTIndices.bin"; SphereInterpolator::SphereInterpolator(const std::string& lutPath) { m_ValidState = true; if (lutPath.length()==0) { if (!LoadLookuptables()) { m_ValidState = false; return; } } else { if (!LoadLookuptables(lutPath)) { m_ValidState = false; return; } } size = 301; sN = (size-1)/2; nverts = ODF_SAMPLING_SIZE; beta = 0.5; inva = (sqrt(1+beta)-sqrt(beta)); b = 1/(1-sqrt(1/beta + 1)); } SphereInterpolator::~SphereInterpolator() { } bool SphereInterpolator::LoadLookuptables(const std::string& lutPath) { MITK_INFO << "SphereInterpolator: loading lookuptables from custom path: " << lutPath; std::string path = lutPath; path.append(BaryCoordsFileName); std::ifstream BaryCoordsStream; BaryCoordsStream.open(path.c_str(), ios::in | ios::binary); MITK_INFO << "SphereInterpolator: 1 " << path; if (!BaryCoordsStream.is_open()) { MITK_INFO << "SphereInterpolator: could not load FiberTrackingLUTBaryCoords.bin from " << path; return false; } - ifstream IndicesStream; + std::ifstream IndicesStream; path = lutPath; path.append("FiberTrackingLUTIndices.bin"); IndicesStream.open(path.c_str(), ios::in | ios::binary); MITK_INFO << "SphereInterpolator: 1 " << path; if (!IndicesStream.is_open()) { MITK_INFO << "SphereInterpolator: could not load FiberTrackingLUTIndices.bin from " << path; return false; } if (LoadLookuptables(BaryCoordsStream, IndicesStream)) { MITK_INFO << "SphereInterpolator: first and second lut loaded successfully"; return true; } return false; } bool SphereInterpolator::LoadLookuptables() { MITK_INFO << "SphereInterpolator: loading lookuptables"; us::Module* module = us::GetModuleContext()->GetModule(); us::ModuleResource BaryCoordsRes = module->GetResource(BaryCoordsFileName); if (!BaryCoordsRes.IsValid()) { MITK_INFO << "Could not retrieve resource " << BaryCoordsFileName; return false; } us::ModuleResource IndicesRes = module->GetResource(IndicesFileName); if (!IndicesRes) { MITK_INFO << "Could not retrieve resource " << IndicesFileName; return false; } us::ModuleResourceStream BaryCoordsStream(BaryCoordsRes, std::ios_base::binary); us::ModuleResourceStream IndicesStream(IndicesRes, std::ios_base::binary); return LoadLookuptables(BaryCoordsStream, IndicesStream); } bool SphereInterpolator::LoadLookuptables(std::istream& BaryCoordsStream, std::istream& IndicesStream) { if (BaryCoordsStream) { try { float tmp; BaryCoordsStream.seekg (0, ios::beg); while (!BaryCoordsStream.eof()) { BaryCoordsStream.read((char *)&tmp, sizeof(tmp)); barycoords.push_back(tmp); } } catch (const std::exception& e) { MITK_INFO << e.what(); } } else { MITK_INFO << "SphereInterpolator: could not load FiberTrackingLUTBaryCoords.bin"; return false; } if (IndicesStream) { try { int tmp; IndicesStream.seekg (0, ios::beg); while (!IndicesStream.eof()) { IndicesStream.read((char *)&tmp, sizeof(tmp)); indices.push_back(tmp); } } catch (const std::exception& e) { MITK_INFO << e.what(); } } else { MITK_INFO << "SphereInterpolator: could not load FiberTrackingLUTIndices.bin"; return false; } return true; } diff --git a/Modules/FiberTracking/Algorithms/itkGibbsTrackingFilter.cpp b/Modules/FiberTracking/Algorithms/itkGibbsTrackingFilter.cpp index fb642db..17d979c 100644 --- a/Modules/FiberTracking/Algorithms/itkGibbsTrackingFilter.cpp +++ b/Modules/FiberTracking/Algorithms/itkGibbsTrackingFilter.cpp @@ -1,515 +1,511 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "itkGibbsTrackingFilter.h" // MITK #include #include #include #include #include //#include #include #include // ITK #include #include #include // MISC #include // #include -#include +#include #include #include #include namespace itk{ template< class ItkOdfImageType > GibbsTrackingFilter< ItkOdfImageType >::GibbsTrackingFilter(): m_StartTemperature(0.1), m_EndTemperature(0.001), m_Iterations(1e9), m_CurrentIteration(0.0), m_CurrentStep(0), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), m_ConnectionPotential(10), m_InexBalance(0), m_ParticlePotential(0.2), m_MinFiberLength(10), m_AbortTracking(false), m_NumAcceptedFibers(0), m_BuildFibers(false), m_ProposalAcceptance(0), m_CurvatureThreshold(0.7), m_DuplicateImage(true), m_NumParticles(0), m_NumConnections(0), m_RandomSeed(-1), m_LoadParameterFile(""), m_LutPath(""), m_IsInValidState(true) { } template< class ItkOdfImageType > GibbsTrackingFilter< ItkOdfImageType >::~GibbsTrackingFilter() { } // fill output fiber bundle datastructure template< class ItkOdfImageType > typename GibbsTrackingFilter< ItkOdfImageType >::FiberPolyDataType GibbsTrackingFilter< ItkOdfImageType >::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } return m_FiberPolyData; } template< class ItkOdfImageType > void GibbsTrackingFilter< ItkOdfImageType > ::EstimateParticleWeight() { MITK_INFO << "GibbsTrackingFilter: estimating particle weight"; float minSpacing; if(m_OdfImage->GetSpacing()[0]GetSpacing()[1] && m_OdfImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_OdfImage->GetSpacing()[0]; else if (m_OdfImage->GetSpacing()[1] < m_OdfImage->GetSpacing()[2]) minSpacing = m_OdfImage->GetSpacing()[1]; else minSpacing = m_OdfImage->GetSpacing()[2]; float m_ParticleLength = 1.5*minSpacing; float m_ParticleWidth = 0.5*minSpacing; // seed random generators Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = Statistics::MersenneTwisterRandomVariateGenerator::New(); if (m_RandomSeed>-1) randGen->SetSeed(m_RandomSeed); else randGen->SetSeed(); // instantiate all necessary components SphereInterpolator* interpolator = new SphereInterpolator(m_LutPath); // handle lookup table not found cases if( !interpolator->IsInValidState() ) { m_IsInValidState = false; m_AbortTracking = true; m_BuildFibers = false; mitkThrow() << "Unable to load lookup tables."; } ParticleGrid* particleGrid = new ParticleGrid(m_MaskImage, m_ParticleLength, m_ParticleGridCellCapacity); GibbsEnergyComputer* encomp = new GibbsEnergyComputer(m_OdfImage, m_MaskImage, particleGrid, interpolator, randGen); // EnergyComputer* encomp = new EnergyComputer(m_OdfImage, m_MaskImage, particleGrid, interpolator, randGen); MetropolisHastingsSampler* sampler = new MetropolisHastingsSampler(particleGrid, encomp, randGen, m_CurvatureThreshold); float alpha = log(m_EndTemperature/m_StartTemperature); m_ParticleWeight = 0.01; int ppv = 0; // main loop int neededParts = 3000; while (ppvSetParameters(m_ParticleWeight,m_ParticleWidth,m_ConnectionPotential*m_ParticleLength*m_ParticleLength,m_CurvatureThreshold,m_InexBalance,m_ParticlePotential); for( int step = 0; step < 10; step++ ) { // update temperatur for simulated annealing process float temperature = m_StartTemperature * exp(alpha*(((1.0)*step)/((1.0)*10))); sampler->SetTemperature(temperature); for (unsigned long i=0; i<10000; i++) sampler->MakeProposal(); } ppv = particleGrid->m_NumParticles; particleGrid->ResetGrid(); } delete sampler; delete encomp; delete particleGrid; delete interpolator; MITK_INFO << "GibbsTrackingFilter: finished estimating particle weight"; } // perform global tracking template< class ItkOdfImageType > void GibbsTrackingFilter< ItkOdfImageType >::GenerateData() { TimeProbe preClock; preClock.Start(); // check if input is Odf or tensor image and generate Odf if necessary if (m_OdfImage.IsNull() && m_TensorImage.IsNotNull()) { TensorImageToOdfImageFilter::Pointer filter = TensorImageToOdfImageFilter::New(); filter->SetInput( m_TensorImage ); filter->Update(); m_OdfImage = filter->GetOutput(); } else if (m_DuplicateImage) // generate local working copy of Odf image (if not disabled) { typedef itk::ImageDuplicator< ItkOdfImageType > DuplicateFilterType; typename DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); duplicator->SetInputImage( m_OdfImage ); duplicator->Update(); m_OdfImage = duplicator->GetOutput(); } // perform mean subtraction on odfs typedef ImageRegionIterator< ItkOdfImageType > InputIteratorType; InputIteratorType it(m_OdfImage, m_OdfImage->GetLargestPossibleRegion() ); it.GoToBegin(); while (!it.IsAtEnd()) { itk::OrientationDistributionFunction odf(it.Get().GetDataPointer()); float mean = odf.GetMeanValue(); odf -= mean; it.Set(odf.GetDataPointer()); ++it; } // check if mask image is given if it needs resampling PrepareMaskImage(); // load parameter file LoadParameters(); // prepare parameters float minSpacing; if(m_OdfImage->GetSpacing()[0]GetSpacing()[1] && m_OdfImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_OdfImage->GetSpacing()[0]; else if (m_OdfImage->GetSpacing()[1] < m_OdfImage->GetSpacing()[2]) minSpacing = m_OdfImage->GetSpacing()[1]; else minSpacing = m_OdfImage->GetSpacing()[2]; if(m_ParticleLength == 0) m_ParticleLength = 1.5*minSpacing; if(m_ParticleWidth == 0) m_ParticleWidth = 0.5*minSpacing; if(m_ParticleWeight == 0) EstimateParticleWeight(); float alpha = log(m_EndTemperature/m_StartTemperature); if (m_CurvatureThreshold < mitk::eps) m_CurvatureThreshold = 0; // seed random generators Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = Statistics::MersenneTwisterRandomVariateGenerator::New(); if (m_RandomSeed>-1) randGen->SetSeed(m_RandomSeed); else randGen->SetSeed(); // load sphere interpolator to evaluate the ODFs SphereInterpolator* interpolator = new SphereInterpolator(m_LutPath); // handle lookup table not found cases if( !interpolator->IsInValidState() ) { m_IsInValidState = false; m_AbortTracking = true; m_BuildFibers = false; mitkThrow() << "Unable to load lookup tables."; } // initialize the actual tracking components (ParticleGrid, Metropolis Hastings Sampler and Energy Computer) ParticleGrid* particleGrid; GibbsEnergyComputer* encomp; MetropolisHastingsSampler* sampler; try{ particleGrid = new ParticleGrid(m_MaskImage, m_ParticleLength, m_ParticleGridCellCapacity); encomp = new GibbsEnergyComputer(m_OdfImage, m_MaskImage, particleGrid, interpolator, randGen); encomp->SetParameters(m_ParticleWeight,m_ParticleWidth,m_ConnectionPotential*m_ParticleLength*m_ParticleLength,m_CurvatureThreshold,m_InexBalance,m_ParticlePotential); sampler = new MetropolisHastingsSampler(particleGrid, encomp, randGen, m_CurvatureThreshold); } catch(...) { MITK_ERROR << "Particle grid allocation failed. Not enough memory? Try to increase the particle length."; m_IsInValidState = false; m_AbortTracking = true; m_BuildFibers = false; return; } MITK_INFO << "----------------------------------------"; MITK_INFO << "Iterations: " << m_Iterations; MITK_INFO << "Particle length: " << m_ParticleLength; MITK_INFO << "Particle width: " << m_ParticleWidth; MITK_INFO << "Particle weight: " << m_ParticleWeight; MITK_INFO << "Start temperature: " << m_StartTemperature; MITK_INFO << "End temperature: " << m_EndTemperature; MITK_INFO << "In/Ex balance: " << m_InexBalance; MITK_INFO << "Min. fiber length: " << m_MinFiberLength; MITK_INFO << "Curvature threshold: " << m_CurvatureThreshold; MITK_INFO << "Random seed: " << m_RandomSeed; MITK_INFO << "----------------------------------------"; // main loop preClock.Stop(); TimeProbe clock; clock.Start(); m_NumAcceptedFibers = 0; m_CurrentIteration = 0; bool just_built_fibers = false; boost::progress_display disp(m_Iterations); if (!m_AbortTracking) while (m_CurrentIterationSetTemperature(temperature); sampler->MakeProposal(); m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/m_CurrentIteration; m_NumParticles = particleGrid->m_NumParticles; m_NumConnections = particleGrid->m_NumConnections; if (m_AbortTracking) break; if (m_BuildFibers) { FiberBuilder fiberBuilder(particleGrid, m_MaskImage); m_FiberPolyData = fiberBuilder.iterate(m_MinFiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); m_BuildFibers = false; just_built_fibers = true; } } if (!just_built_fibers) { FiberBuilder fiberBuilder(particleGrid, m_MaskImage); m_FiberPolyData = fiberBuilder.iterate(m_MinFiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); } clock.Stop(); delete sampler; delete encomp; delete interpolator; delete particleGrid; m_AbortTracking = true; m_BuildFibers = false; int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "GibbsTrackingFilter: finished gibbs tracking in " << h << "h, " << m << "m and " << s << "s"; m = (int)preClock.GetTotal()/60; s = (int)preClock.GetTotal()%60; MITK_INFO << "GibbsTrackingFilter: preparation of the data took " << m << "m and " << s << "s"; MITK_INFO << "GibbsTrackingFilter: " << m_NumAcceptedFibers << " fibers accepted"; // sampler->PrintProposalTimes(); SaveParameters(); } template< class ItkOdfImageType > void GibbsTrackingFilter< ItkOdfImageType >::PrepareMaskImage() { if(m_MaskImage.IsNull()) { MITK_INFO << "GibbsTrackingFilter: generating default mask image"; m_MaskImage = ItkFloatImageType::New(); m_MaskImage->SetSpacing( m_OdfImage->GetSpacing() ); m_MaskImage->SetOrigin( m_OdfImage->GetOrigin() ); m_MaskImage->SetDirection( m_OdfImage->GetDirection() ); m_MaskImage->SetRegions( m_OdfImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1.0); } else if ( m_MaskImage->GetLargestPossibleRegion().GetSize()[0]!=m_OdfImage->GetLargestPossibleRegion().GetSize()[0] || m_MaskImage->GetLargestPossibleRegion().GetSize()[1]!=m_OdfImage->GetLargestPossibleRegion().GetSize()[1] || m_MaskImage->GetLargestPossibleRegion().GetSize()[2]!=m_OdfImage->GetLargestPossibleRegion().GetSize()[2] || m_MaskImage->GetSpacing()[0]!=m_OdfImage->GetSpacing()[0] || m_MaskImage->GetSpacing()[1]!=m_OdfImage->GetSpacing()[1] || m_MaskImage->GetSpacing()[2]!=m_OdfImage->GetSpacing()[2] ) { MITK_INFO << "GibbsTrackingFilter: resampling mask image"; typedef itk::ResampleImageFilter< ItkFloatImageType, ItkFloatImageType, float > ResamplerType; ResamplerType::Pointer resampler = ResamplerType::New(); resampler->SetOutputSpacing( m_OdfImage->GetSpacing() ); resampler->SetOutputOrigin( m_OdfImage->GetOrigin() ); resampler->SetOutputDirection( m_OdfImage->GetDirection() ); resampler->SetSize( m_OdfImage->GetLargestPossibleRegion().GetSize() ); resampler->SetInput( m_MaskImage ); resampler->SetDefaultPixelValue(0.0); resampler->Update(); m_MaskImage = resampler->GetOutput(); MITK_INFO << "GibbsTrackingFilter: resampling finished"; } } // load tracking paramters from xml file (.gtp) template< class ItkOdfImageType > bool GibbsTrackingFilter< ItkOdfImageType >::LoadParameters() { m_AbortTracking = true; try { if( m_LoadParameterFile.length()==0 ) { m_AbortTracking = false; return true; } MITK_INFO << "GibbsTrackingFilter: loading parameter file " << m_LoadParameterFile; - TiXmlDocument doc( m_LoadParameterFile ); - doc.LoadFile(); + tinyxml2::XMLDocument doc; + doc.LoadFile(m_LoadParameterFile.c_str()); - TiXmlHandle hDoc(&doc); - TiXmlElement* pElem; - TiXmlHandle hRoot(nullptr); - - pElem = hDoc.FirstChildElement().Element(); - hRoot = TiXmlHandle(pElem); - pElem = hRoot.FirstChildElement("parameter_set").Element(); + + tinyxml2::XMLNode * hRoot = doc.FirstChild(); + + auto pElem = hRoot->FirstChildElement("parameter_set"); std::string iterations(pElem->Attribute("iterations")); m_Iterations = boost::lexical_cast(iterations); std::string particleLength(pElem->Attribute("particle_length")); m_ParticleLength = boost::lexical_cast(particleLength); std::string particleWidth(pElem->Attribute("particle_width")); m_ParticleWidth = boost::lexical_cast(particleWidth); std::string partWeight(pElem->Attribute("particle_weight")); m_ParticleWeight = boost::lexical_cast(partWeight); std::string startTemp(pElem->Attribute("temp_start")); m_StartTemperature = boost::lexical_cast(startTemp); std::string endTemp(pElem->Attribute("temp_end")); m_EndTemperature = boost::lexical_cast(endTemp); std::string inExBalance(pElem->Attribute("inexbalance")); m_InexBalance = boost::lexical_cast(inExBalance); std::string fiberLength(pElem->Attribute("fiber_length")); m_MinFiberLength = boost::lexical_cast(fiberLength); std::string curvThres(pElem->Attribute("curvature_threshold")); m_CurvatureThreshold = cos(boost::lexical_cast(curvThres)*itk::Math::pi/180); m_AbortTracking = false; MITK_INFO << "GibbsTrackingFilter: parameter file loaded successfully"; return true; } catch(...) { MITK_INFO << "GibbsTrackingFilter: could not load parameter file"; return false; } } // save current tracking paramters to xml file (.gtp) template< class ItkOdfImageType > bool GibbsTrackingFilter< ItkOdfImageType >::SaveParameters() { try { if( m_SaveParameterFile.length()==0 ) { MITK_INFO << "GibbsTrackingFilter: no filename specified to save parameters"; return true; } MITK_INFO << "GibbsTrackingFilter: saving parameter file " << m_SaveParameterFile; - - TiXmlDocument documentXML; - TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); - documentXML.LinkEndChild( declXML ); - - TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); - mainXML->SetAttribute("file_version", "0.1"); - documentXML.LinkEndChild(mainXML); - - TiXmlElement* paramXML = new TiXmlElement("parameter_set"); - paramXML->SetAttribute("iterations", boost::lexical_cast(m_Iterations)); - paramXML->SetAttribute("particle_length", boost::lexical_cast(m_ParticleLength)); - paramXML->SetAttribute("particle_width", boost::lexical_cast(m_ParticleWidth)); - paramXML->SetAttribute("particle_weight", boost::lexical_cast(m_ParticleWeight)); - paramXML->SetAttribute("temp_start", boost::lexical_cast(m_StartTemperature)); - paramXML->SetAttribute("temp_end", boost::lexical_cast(m_EndTemperature)); - paramXML->SetAttribute("inexbalance", boost::lexical_cast(m_InexBalance)); - paramXML->SetAttribute("fiber_length", boost::lexical_cast(m_MinFiberLength)); - paramXML->SetAttribute("curvature_threshold", boost::lexical_cast(m_CurvatureThreshold)); - mainXML->LinkEndChild(paramXML); + + tinyxml2::XMLDocument documentXML; + //TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); + //documentXML.LinkEndChild( declXML ); + + tinyxml2::XMLNode * mainXML = documentXML.NewElement("global_tracking_parameter_file"); + documentXML.InsertFirstChild(mainXML); + + auto paramXML = documentXML.NewElement("parameter_set"); + paramXML->SetAttribute("iterations", m_Iterations); + paramXML->SetAttribute("particle_length", m_ParticleLength); + paramXML->SetAttribute("particle_width", m_ParticleWidth); + paramXML->SetAttribute("particle_weight", m_ParticleWeight); + paramXML->SetAttribute("temp_start", m_StartTemperature); + paramXML->SetAttribute("temp_end", m_EndTemperature); + paramXML->SetAttribute("inexbalance", m_InexBalance); + paramXML->SetAttribute("fiber_length", m_MinFiberLength); + paramXML->SetAttribute("curvature_threshold", m_CurvatureThreshold); + mainXML->InsertEndChild(paramXML); if(!boost::algorithm::ends_with(m_SaveParameterFile, ".gtp")) m_SaveParameterFile.append(".gtp"); - documentXML.SaveFile( m_SaveParameterFile ); + documentXML.SaveFile( m_SaveParameterFile.c_str() ); MITK_INFO << "GibbsTrackingFilter: parameter file saved successfully"; return true; } catch(...) { MITK_INFO << "GibbsTrackingFilter: could not save parameter file"; return false; } } template< class ItkOdfImageType > void GibbsTrackingFilter< ItkOdfImageType >::SetDicomProperties(mitk::FiberBundle::Pointer fib) { std::string model_code_value = "-"; std::string model_code_meaning = "-"; std::string algo_code_value = "sup181_ee03"; std::string algo_code_meaning = "Global"; fib->SetProperty("DICOM.anatomy.value", mitk::StringProperty::New("T-A0095")); fib->SetProperty("DICOM.anatomy.meaning", mitk::StringProperty::New("White matter of brain and spinal cord")); fib->SetProperty("DICOM.algo_family_code.value", mitk::StringProperty::New(algo_code_value)); fib->SetProperty("DICOM.algo_family_code.meaning", mitk::StringProperty::New(algo_code_meaning)); fib->SetProperty("DICOM.model_code.value", mitk::StringProperty::New(model_code_value)); fib->SetProperty("DICOM.model_code.meaning", mitk::StringProperty::New(model_code_meaning)); } } diff --git a/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.h b/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.h index 5607ba4..da734ad 100644 --- a/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.h +++ b/Modules/MriSimulation/Algorithms/itkTractsToDWIImageFilter.h @@ -1,177 +1,177 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkTractsToDWIImageFilter_h__ #define __itkTractsToDWIImageFilter_h__ #include #include #include #include #include #include #include #include #include #include namespace itk { /** * \brief Generates artificial diffusion weighted image volume from the input fiberbundle using a generic multicompartment model. * See "Fiberfox: Facilitating the creation of realistic white matter software phantoms" (DOI: 10.1002/mrm.25045) for details. */ template< class PixelType > class TractsToDWIImageFilter : public ImageSource< itk::VectorImage< PixelType, 3 > > { public: typedef TractsToDWIImageFilter Self; typedef ImageSource< itk::VectorImage< PixelType, 3 > > Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef typename Superclass::OutputImageType OutputImageType; typedef itk::Image ItkDoubleImgType4D; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image ItkUcharImgType; typedef mitk::FiberBundle::Pointer FiberBundleType; typedef itk::VectorImage< double, 3 > DoubleDwiType; typedef itk::Matrix MatrixType; typedef itk::Image< float, 2 > Float2DImageType; typedef itk::Image< vcl_complex< float >, 2 > Complex2DImageType; typedef itk::Vector< float, 3> VectorType; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( TractsToDWIImageFilter, ImageSource ) /** Input */ itkSetMacro( FiberBundle, FiberBundleType ) ///< Input fiber bundle itkSetMacro( InputImage, typename OutputImageType::Pointer ) ///< Input diffusion-weighted image. If no fiber bundle is set, then the acquisition is simulated for this image without a new diffusion simulation. itkSetMacro( UseConstantRandSeed, bool ) ///< Seed for random generator. void SetParameters( FiberfoxParameters param ) ///< Simulation parameters. { m_Parameters = param; } /** Output */ FiberfoxParameters GetParameters(){ return m_Parameters; } std::vector< ItkDoubleImgType::Pointer > GetVolumeFractions() ///< one double image for each compartment containing the corresponding volume fraction per voxel { return m_VolumeFractions; } itkGetMacro( StatusText, std::string ) itkGetMacro( PhaseImage, DoubleDwiType::Pointer ) itkGetMacro( KspaceImage, DoubleDwiType::Pointer ) itkGetMacro( CoilPointset, mitk::PointSet::Pointer ) itkGetMacro( TickImage, typename Float2DImageType::Pointer ) ///< k-space readout ordering encoded in the voxels itkGetMacro( RfImage, typename Float2DImageType::Pointer ) ///< time passed since last RF pulse encoded per voxel void GenerateData() override; std::vector GetOutputImagesReal() const { return m_OutputImagesReal; } std::vector GetOutputImagesImag() const { return m_OutputImagesImag; } protected: TractsToDWIImageFilter(); ~TractsToDWIImageFilter() override; itk::Vector GetItkVector(double point[3]); vnl_vector_fixed GetVnlVector(double point[3]); vnl_vector_fixed GetVnlVector(Vector< float, 3 >& vector); double RoundToNearest(double num); std::string GetTime(); bool PrepareLogFile(); /** Prepares the log file and returns true if successful or false if failed. */ void PrintToLog(std::string m, bool addTime=true, bool linebreak=true, bool stdOut=true); /** Transform generated image compartment by compartment, channel by channel and slice by slice using DFT and add k-space artifacts/effects. */ DoubleDwiType::Pointer SimulateKspaceAcquisition(std::vector< DoubleDwiType::Pointer >& images); /** Generate signal of non-fiber compartments. */ void SimulateExtraAxonalSignal(ItkUcharImgType::IndexType& index, itk::Point& volume_fraction_point, double intraAxonalVolume, int g); /** Move fibers to simulate headmotion */ void SimulateMotion(int g=-1); void CheckVolumeFractionImages(); ItkDoubleImgType::Pointer NormalizeInsideMask(ItkDoubleImgType::Pointer image); void InitializeData(); void InitializeFiberData(); itk::Point GetMovedPoint(itk::Index<3>& index, bool forward); // input mitk::FiberfoxParameters m_Parameters; FiberBundleType m_FiberBundle; typename OutputImageType::Pointer m_InputImage; // output typename OutputImageType::Pointer m_OutputImage; typename DoubleDwiType::Pointer m_PhaseImage; typename DoubleDwiType::Pointer m_KspaceImage; std::vector < typename DoubleDwiType::Pointer > m_OutputImagesReal; std::vector < typename DoubleDwiType::Pointer > m_OutputImagesImag; std::vector< ItkDoubleImgType::Pointer > m_VolumeFractions; std::string m_StatusText; // MISC itk::TimeProbe m_TimeProbe; bool m_UseConstantRandSeed; bool m_MaskImageSet; - ofstream m_Logfile; + std::ofstream m_Logfile; std::string m_MotionLog; std::string m_SpikeLog; // signal generation FiberBundleType m_FiberBundleTransformed; ///< transformed bundle simulating headmotion itk::Vector m_WorkingSpacing; itk::Point m_WorkingOrigin; ImageRegion<3> m_WorkingImageRegion; double m_VoxelVolume; std::vector< DoubleDwiType::Pointer > m_CompartmentImages; ItkUcharImgType::Pointer m_TransformedMaskImage; ///< copy of mask image (changes for each motion step) ItkUcharImgType::Pointer m_UpsampledMaskImage; ///< helper image for motion simulation std::vector< MatrixType > m_RotationsInv; std::vector< MatrixType > m_Rotations; /// m_Translations; ///::Pointer m_DoubleInterpolator; itk::Vector m_NullDir; Float2DImageType::Pointer m_TickImage; Float2DImageType::Pointer m_RfImage; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkTractsToDWIImageFilter.cpp" #endif #endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/src/internal/QmitkPartialVolumeAnalysisView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/src/internal/QmitkPartialVolumeAnalysisView.cpp index c808f0d..88e6e71 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/src/internal/QmitkPartialVolumeAnalysisView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.partialvolume/src/internal/QmitkPartialVolumeAnalysisView.cpp @@ -1,2139 +1,2139 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkPartialVolumeAnalysisView.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include "QmitkSliderNavigatorWidget.h" #include #include "mitkNodePredicateDataType.h" #include "mitkNodePredicateOr.h" #include "mitkImageTimeSelector.h" #include "mitkProperties.h" #include "mitkProgressBar.h" #include "mitkImageCast.h" #include "mitkImageToItk.h" #include "mitkITKImageImport.h" #include "mitkDataNodeObject.h" #include "mitkNodePredicateData.h" #include "mitkPlanarFigureInteractor.h" #include "mitkTensorImage.h" #include "mitkPlanarCircle.h" #include "mitkPlanarRectangle.h" #include "mitkPlanarPolygon.h" #include "mitkPartialVolumeAnalysisClusteringCalculator.h" #include "usModuleRegistry.h" #include #include "itkTensorDerivedMeasurementsFilter.h" #include "itkDiffusionTensor3D.h" #include "itkCartesianToPolarVectorImageFilter.h" #include "itkPolarToCartesianVectorImageFilter.h" #include "itkBinaryThresholdImageFilter.h" #include "itkMaskImageFilter.h" #include "itkCastImageFilter.h" #include "itkImageMomentsCalculator.h" #include #include #include #include #include const std::string QmitkPartialVolumeAnalysisView::VIEW_ID = "org.mitk.views.partialvolumeanalysisview"; class QmitkRequestStatisticsUpdateEvent : public QEvent { public: enum Type { StatisticsUpdateRequest = QEvent::MaxUser - 1025 }; QmitkRequestStatisticsUpdateEvent() : QEvent( (QEvent::Type) StatisticsUpdateRequest ) {}; }; typedef itk::Image ImageType; typedef itk::Image FloatImageType; typedef itk::Image, 3> VectorImageType; inline bool my_isnan(float x) { volatile float d = x; if(d!=d) return true; if(d==d) return false; return d != d; } QmitkPartialVolumeAnalysisView::QmitkPartialVolumeAnalysisView(QObject * /*parent*/, const char * /*name*/) : QmitkAbstractView(), m_Controls(nullptr), m_TimeStepperAdapter(nullptr), m_MeasurementInfoRenderer(nullptr), m_MeasurementInfoAnnotation(nullptr), m_IsTensorImage(false), m_SelectedRenderWindow(nullptr), m_LastRenderWindow(nullptr), m_ImageObserverTag(-1), m_ImageMaskObserverTag(-1), m_PlanarFigureObserverTag(-1), m_CurrentStatisticsValid(false), m_StatisticsUpdatePending(false), m_GaussianSigmaChangedSliding(false), m_NumberBinsSliding(false), m_UpsamplingChangedSliding(false), m_EllipseCounter(0), m_RectangleCounter(0), m_PolygonCounter(0), m_CurrentFigureNodeInitialized(false), m_QuantifyClass(2), m_IconTexOFF(new QIcon(":/QmitkPartialVolumeAnalysisView/texIntOFFIcon.png")), m_IconTexON(new QIcon(":/QmitkPartialVolumeAnalysisView/texIntONIcon.png")), m_TexIsOn(true), m_Visible(false) { } QmitkPartialVolumeAnalysisView::~QmitkPartialVolumeAnalysisView() { if ( m_SelectedImage.IsNotNull() ) m_SelectedImage->RemoveObserver( m_ImageObserverTag ); if ( m_SelectedImageMask.IsNotNull() ) m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); if ( m_SelectedPlanarFigure.IsNotNull() ) { m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); m_SelectedPlanarFigure->RemoveObserver( m_InitializedObserverTag ); } this->GetDataStorage()->AddNodeEvent -= mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeAddedInDataStorage ); m_SelectedPlanarFigureNodes->NodeChanged.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); m_SelectedPlanarFigureNodes->NodeRemoved.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); m_SelectedPlanarFigureNodes->PropertyChanged.RemoveListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); m_SelectedImageNodes->NodeChanged.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); m_SelectedImageNodes->NodeRemoved.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); m_SelectedImageNodes->PropertyChanged.RemoveListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); } void QmitkPartialVolumeAnalysisView::CreateQtPartControl(QWidget *parent) { if (m_Controls == nullptr) { m_Controls = new Ui::QmitkPartialVolumeAnalysisViewControls; m_Controls->setupUi(parent); this->CreateConnections(); } SetHistogramVisibility(); m_Controls->m_TextureIntON->setIcon(*m_IconTexON); m_Controls->m_SimilarAnglesFrame->setVisible(false); m_Controls->m_SimilarAnglesLabel->setVisible(false); vtkTextProperty *textProp = vtkTextProperty::New(); textProp->SetColor(1.0, 1.0, 1.0); m_MeasurementInfoAnnotation = vtkCornerAnnotation::New(); m_MeasurementInfoAnnotation->SetMaximumFontSize(12); m_MeasurementInfoAnnotation->SetTextProperty(textProp); m_MeasurementInfoRenderer = vtkRenderer::New(); m_MeasurementInfoRenderer->AddActor(m_MeasurementInfoAnnotation); m_SelectedPlanarFigureNodes = mitk::DataStorageSelection::New(this->GetDataStorage(), false); m_SelectedPlanarFigureNodes->NodeChanged.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); m_SelectedPlanarFigureNodes->NodeRemoved.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); m_SelectedPlanarFigureNodes->PropertyChanged.AddListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); m_SelectedImageNodes = mitk::DataStorageSelection::New(this->GetDataStorage(), false); m_SelectedImageNodes->PropertyChanged.AddListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); m_SelectedImageNodes->NodeChanged.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); m_SelectedImageNodes->NodeRemoved.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); this->GetDataStorage()->AddNodeEvent.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeAddedInDataStorage ) ); Select(nullptr,true,true); SetAdvancedVisibility(); } void QmitkPartialVolumeAnalysisView::SetFocus() { m_Controls->m_CircleButton->setFocus(); } void QmitkPartialVolumeAnalysisView::SetHistogramVisibility() { m_Controls->m_HistogramWidget->setVisible(m_Controls->m_DisplayHistogramCheckbox->isChecked()); } void QmitkPartialVolumeAnalysisView::SetAdvancedVisibility() { m_Controls->frame_7->setVisible(m_Controls->m_AdvancedCheckbox->isChecked()); } void QmitkPartialVolumeAnalysisView::CreateConnections() { if ( m_Controls ) { connect( m_Controls->m_DisplayHistogramCheckbox, SIGNAL( clicked() ) , this, SLOT( SetHistogramVisibility() ) ); connect( m_Controls->m_AdvancedCheckbox, SIGNAL( clicked() ) , this, SLOT( SetAdvancedVisibility() ) ); connect( m_Controls->m_NumberBinsSlider, SIGNAL( sliderReleased () ), this, SLOT( NumberBinsReleasedSlider( ) ) ); connect( m_Controls->m_UpsamplingSlider, SIGNAL( sliderReleased( ) ), this, SLOT( UpsamplingReleasedSlider( ) ) ); connect( m_Controls->m_GaussianSigmaSlider, SIGNAL( sliderReleased( ) ), this, SLOT( GaussianSigmaReleasedSlider( ) ) ); connect( m_Controls->m_SimilarAnglesSlider, SIGNAL( sliderReleased( ) ), this, SLOT( SimilarAnglesReleasedSlider( ) ) ); connect( m_Controls->m_NumberBinsSlider, SIGNAL( valueChanged (int) ), this, SLOT( NumberBinsChangedSlider( int ) ) ); connect( m_Controls->m_UpsamplingSlider, SIGNAL( valueChanged( int ) ), this, SLOT( UpsamplingChangedSlider( int ) ) ); connect( m_Controls->m_GaussianSigmaSlider, SIGNAL( valueChanged( int ) ), this, SLOT( GaussianSigmaChangedSlider( int ) ) ); connect( m_Controls->m_SimilarAnglesSlider, SIGNAL( valueChanged( int ) ), this, SLOT( SimilarAnglesChangedSlider(int) ) ); connect( m_Controls->m_OpacitySlider, SIGNAL( valueChanged( int ) ), this, SLOT( OpacityChangedSlider(int) ) ); connect( (QObject*)(m_Controls->m_ButtonCopyHistogramToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(ToClipBoard())); connect( m_Controls->m_CircleButton, SIGNAL( clicked() ) , this, SLOT( ActionDrawEllipseTriggered() ) ); connect( m_Controls->m_RectangleButton, SIGNAL( clicked() ) , this, SLOT( ActionDrawRectangleTriggered() ) ); connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ) , this, SLOT( ActionDrawPolygonTriggered() ) ); connect( m_Controls->m_GreenRadio, SIGNAL( clicked(bool) ) , this, SLOT( GreenRadio(bool) ) ); connect( m_Controls->m_PartialVolumeRadio, SIGNAL( clicked(bool) ) , this, SLOT( PartialVolumeRadio(bool) ) ); connect( m_Controls->m_BlueRadio, SIGNAL( clicked(bool) ) , this, SLOT( BlueRadio(bool) ) ); connect( m_Controls->m_AllRadio, SIGNAL( clicked(bool) ) , this, SLOT( AllRadio(bool) ) ); connect( m_Controls->m_EstimateCircle, SIGNAL( clicked() ) , this, SLOT( EstimateCircle() ) ); connect( (QObject*)(m_Controls->m_TextureIntON), SIGNAL(clicked()), this, SLOT(TextIntON()) ); connect( m_Controls->m_ExportClusteringResultsButton, SIGNAL(clicked()), this, SLOT(ExportClusteringResults())); } } void QmitkPartialVolumeAnalysisView::ExportClusteringResults() { if (m_ClusteringResult.IsNull() || m_SelectedImage.IsNull()) return; mitk::BaseGeometry* geometry = m_SelectedImage->GetGeometry(); itk::Image< short, 3>::Pointer referenceImage = itk::Image< short, 3>::New(); itk::Vector newSpacing = geometry->GetSpacing(); mitk::Point3D newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); itk::Matrix newDirection; itk::ImageRegion<3> imageRegion; for (int i=0; i<3; i++) for (int j=0; j<3; j++) newDirection[j][i] = geometry->GetMatrixColumn(i)[j]/newSpacing[j]; imageRegion.SetSize(0, geometry->GetExtent(0)); imageRegion.SetSize(1, geometry->GetExtent(1)); imageRegion.SetSize(2, geometry->GetExtent(2)); // apply new image parameters referenceImage->SetSpacing( newSpacing ); referenceImage->SetOrigin( newOrigin ); referenceImage->SetDirection( newDirection ); referenceImage->SetRegions( imageRegion ); referenceImage->Allocate(); typedef itk::Image< float, 3 > OutType; mitk::Image::Pointer mitkInImage = dynamic_cast(m_ClusteringResult->GetData()); typedef itk::Image< itk::RGBAPixel, 3 > ItkRgbaImageType; typedef mitk::ImageToItk< ItkRgbaImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkInImage); caster->Update(); ItkRgbaImageType::Pointer itkInImage = caster->GetOutput(); typedef itk::ExtractChannelFromRgbaImageFilter< itk::Image< short, 3>, OutType > ExtractionFilterType; ExtractionFilterType::Pointer filter = ExtractionFilterType::New(); filter->SetInput(itkInImage); filter->SetChannel(ExtractionFilterType::ALPHA); filter->SetReferenceImage(referenceImage); filter->Update(); OutType::Pointer outImg = filter->GetOutput(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); // init data node mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); node->SetName("Clustering Result"); GetDataStorage()->Add(node); } void QmitkPartialVolumeAnalysisView::EstimateCircle() { typedef itk::Image SegImageType; SegImageType::Pointer mask_itk = SegImageType::New(); typedef mitk::ImageToItk CastType; CastType::Pointer caster = CastType::New(); caster->SetInput(m_SelectedImageMask); caster->Update(); typedef itk::ImageMomentsCalculator< SegImageType > MomentsType; MomentsType::Pointer momentsCalc = MomentsType::New(); momentsCalc->SetImage(caster->GetOutput()); momentsCalc->Compute(); MomentsType::VectorType cog = momentsCalc->GetCenterOfGravity(); MomentsType::MatrixType axes = momentsCalc->GetPrincipalAxes(); // moments-coord conversion // third coordinate min oder max? // max-min = extent MomentsType::AffineTransformPointer trafo = momentsCalc->GetPhysicalAxesToPrincipalAxesTransform(); itk::ImageRegionIterator itimage(caster->GetOutput(), caster->GetOutput()->GetLargestPossibleRegion()); itimage.GoToBegin(); double max = -9999999999.0; double min = 9999999999.0; while( !itimage.IsAtEnd() ) { if(itimage.Get()) { ImageType::IndexType index = itimage.GetIndex(); itk::Point point; caster->GetOutput()->TransformIndexToPhysicalPoint(index,point); itk::Point newPoint; newPoint = trafo->TransformPoint(point); if(newPoint[2]max) max = newPoint[2]; } ++itimage; } double extent = max - min; MITK_DEBUG << "EXTENT = " << extent; mitk::Point3D origin; mitk::Vector3D right, bottom, normal; double factor = 1000.0; mitk::FillVector3D(origin, cog[0]-factor*axes[1][0]-factor*axes[2][0], cog[1]-factor*axes[1][1]-factor*axes[2][1], cog[2]-factor*axes[1][2]-factor*axes[2][2]); // mitk::FillVector3D(normal, axis[0][0],axis[0][1],axis[0][2]); mitk::FillVector3D(bottom, 2*factor*axes[1][0], 2*factor*axes[1][1], 2*factor*axes[1][2]); mitk::FillVector3D(right, 2*factor*axes[2][0], 2*factor*axes[2][1], 2*factor*axes[2][2]); mitk::PlaneGeometry::Pointer planegeometry = mitk::PlaneGeometry::New(); planegeometry->InitializeStandardPlane(right.GetVnlVector(), bottom.GetVnlVector()); planegeometry->SetOrigin(origin); double len1 = sqrt(axes[1][0]*axes[1][0] + axes[1][1]*axes[1][1] + axes[1][2]*axes[1][2]); double len2 = sqrt(axes[2][0]*axes[2][0] + axes[2][1]*axes[2][1] + axes[2][2]*axes[2][2]); mitk::Point2D point1; point1[0] = factor*len1; point1[1] = factor*len2; mitk::Point2D point2; point2[0] = factor*len1+extent*.5; point2[1] = factor*len2; mitk::PlanarCircle::Pointer circle = mitk::PlanarCircle::New(); circle->SetPlaneGeometry(planegeometry); circle->PlaceFigure( point1 ); circle->SetControlPoint(0,point1); circle->SetControlPoint(1,point2); //circle->SetCurrentControlPoint( point2 ); mitk::PlanarFigure::PolyLineType polyline = circle->GetPolyLine( 0 ); MITK_DEBUG << "SIZE of planar figure polyline: " << polyline.size(); AddFigureToDataStorage(circle, "Circle"); } bool QmitkPartialVolumeAnalysisView::AssertDrawingIsPossible(bool checked) { if (m_SelectedImageNodes->GetNode().IsNull()) { checked = false; this->HandleException("Please select an image!", dynamic_cast(this->parent()), true); return false; } //this->GetRenderWindowPart(OPEN)->EnableSlicingPlanes(false); return checked; } void QmitkPartialVolumeAnalysisView::ActionDrawEllipseTriggered() { bool checked = m_Controls->m_CircleButton->isChecked(); if(!this->AssertDrawingIsPossible(checked)) return; mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); // using PV_ prefix for planar figures from this view // to distinguish them from that ones created throught the measurement view this->AddFigureToDataStorage(figure, QString("PV_Circle%1").arg(++m_EllipseCounter)); MITK_DEBUG << "PlanarCircle created ..."; } void QmitkPartialVolumeAnalysisView::ActionDrawRectangleTriggered() { bool checked = m_Controls->m_RectangleButton->isChecked(); if(!this->AssertDrawingIsPossible(checked)) return; mitk::PlanarRectangle::Pointer figure = mitk::PlanarRectangle::New(); // using PV_ prefix for planar figures from this view // to distinguish them from that ones created throught the measurement view this->AddFigureToDataStorage(figure, QString("PV_Rectangle%1").arg(++m_RectangleCounter)); MITK_DEBUG << "PlanarRectangle created ..."; } void QmitkPartialVolumeAnalysisView::ActionDrawPolygonTriggered() { bool checked = m_Controls->m_PolygonButton->isChecked(); if(!this->AssertDrawingIsPossible(checked)) return; mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); figure->ClosedOn(); // using PV_ prefix for planar figures from this view // to distinguish them from that ones created throught the measurement view this->AddFigureToDataStorage(figure, QString("PV_Polygon%1").arg(++m_PolygonCounter)); MITK_DEBUG << "PlanarPolygon created ..."; } void QmitkPartialVolumeAnalysisView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, const char *propertyKey, mitk::BaseProperty *property ) { mitk::DataNode::Pointer newNode = mitk::DataNode::New(); newNode->SetName(name.toStdString()); newNode->SetData(figure); // Add custom property, if available if ( (propertyKey != nullptr) && (property != nullptr) ) { newNode->AddProperty( propertyKey, property ); } // figure drawn on the topmost layer / image this->GetDataStorage()->Add(newNode, m_SelectedImageNodes->GetNode() ); QList selectedNodes = this->GetDataManagerSelection(); for(int i = 0; i < selectedNodes.size(); ++i) { selectedNodes[i]->SetSelected(false); } std::vector selectedPFNodes = m_SelectedPlanarFigureNodes->GetNodes(); for(std::size_t i = 0; i < selectedPFNodes.size(); ++i) { selectedPFNodes[i]->SetSelected(false); } newNode->SetSelected(true); Select(newNode); } void QmitkPartialVolumeAnalysisView::PlanarFigureInitialized() { if(m_SelectedPlanarFigureNodes->GetNode().IsNull()) return; m_CurrentFigureNodeInitialized = true; this->Select(m_SelectedPlanarFigureNodes->GetNode()); m_Controls->m_CircleButton->setChecked(false); m_Controls->m_RectangleButton->setChecked(false); m_Controls->m_PolygonButton->setChecked(false); //this->GetRenderWindowPart(OPEN)->EnableSlicingPlanes(true); this->RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::PlanarFigureFocus(mitk::DataNode* node) { mitk::PlanarFigure* _PlanarFigure = 0; _PlanarFigure = dynamic_cast (node->GetData()); if (_PlanarFigure) { FindRenderWindow(node); const mitk::PlaneGeometry* _PlaneGeometry = _PlanarFigure->GetPlaneGeometry(); // make node visible if (m_SelectedRenderWindow) { mitk::Point3D centerP = _PlaneGeometry->GetOrigin(); m_SelectedRenderWindow->GetSliceNavigationController()->ReorientSlices( centerP, _PlaneGeometry->GetNormal()); m_SelectedRenderWindow->GetSliceNavigationController()->SelectSliceByPoint( centerP); } } } void QmitkPartialVolumeAnalysisView::FindRenderWindow(mitk::DataNode* node) { if (node && dynamic_cast (node->GetData())) { m_SelectedRenderWindow = 0; bool PlanarFigureInitializedWindow = false; foreach(QmitkRenderWindow * window, this->GetRenderWindowPart()->GetQmitkRenderWindows().values()) { if (!m_SelectedRenderWindow && node->GetBoolProperty("PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, window->GetRenderer())) { m_SelectedRenderWindow = window; } } } } void QmitkPartialVolumeAnalysisView::OnSelectionChanged(berry::IWorkbenchPart::Pointer, const QList &nodes) { m_Controls->m_InputData->setTitle("Please Select Input Data"); if (!m_Visible) return; if ( nodes.empty() ) { if (m_ClusteringResult.IsNotNull()) { this->GetDataStorage()->Remove(m_ClusteringResult); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } Select(nullptr, true, true); } for (int i=0; iRemoveOrphanImages(); bool somethingChanged = false; if(node.IsNull()) { somethingChanged = true; if(clearMaskOnFirstArgnullptr) { if ( (m_SelectedImageMask.IsNotNull()) && (m_ImageMaskObserverTag >= 0) ) { m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); m_ImageMaskObserverTag = -1; } if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_PlanarFigureObserverTag >= 0) ) { m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); m_PlanarFigureObserverTag = -1; } if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_InitializedObserverTag >= 0) ) { m_SelectedPlanarFigure->RemoveObserver( m_InitializedObserverTag ); m_InitializedObserverTag = -1; } m_SelectedPlanarFigure = nullptr; m_SelectedPlanarFigureNodes->RemoveAllNodes(); m_CurrentFigureNodeInitialized = false; m_SelectedRenderWindow = 0; m_SelectedMaskNode = nullptr; m_SelectedImageMask = nullptr; } if(clearImageOnFirstArgnullptr) { if ( (m_SelectedImage.IsNotNull()) && (m_ImageObserverTag >= 0) ) { m_SelectedImage->RemoveObserver( m_ImageObserverTag ); m_ImageObserverTag = -1; } m_SelectedImageNodes->RemoveAllNodes(); m_SelectedImage = nullptr; m_IsTensorImage = false; m_FAImage = nullptr; m_RDImage = nullptr; m_ADImage = nullptr; m_MDImage = nullptr; m_CAImage = nullptr; m_DirectionComp1Image = nullptr; m_DirectionComp2Image = nullptr; m_AngularErrorImage = nullptr; m_Controls->m_SimilarAnglesFrame->setVisible(false); m_Controls->m_SimilarAnglesLabel->setVisible(false); } } else { typedef itk::SimpleMemberCommand< QmitkPartialVolumeAnalysisView > ITKCommandType; ITKCommandType::Pointer changeListener; changeListener = ITKCommandType::New(); changeListener->SetCallbackFunction( this, &QmitkPartialVolumeAnalysisView::RequestStatisticsUpdate ); // Get selected element mitk::TensorImage *selectedTensorImage = dynamic_cast< mitk::TensorImage * >( node->GetData() ); mitk::Image *selectedImage = dynamic_cast< mitk::Image * >( node->GetData() ); mitk::PlanarFigure *selectedPlanar = dynamic_cast< mitk::PlanarFigure * >( node->GetData() ); bool isMask = false; bool isImage = false; bool isPlanar = false; bool isTensorImage = false; if (selectedTensorImage != nullptr) { isTensorImage = true; } else if(selectedImage != nullptr) { node->GetPropertyValue("binary", isMask); isImage = !isMask; } else if ( (selectedPlanar != nullptr) ) { isPlanar = true; } // image if(isImage && selectedImage->GetDimension()==3) { if(selectedImage != m_SelectedImage.GetPointer()) { somethingChanged = true; if ( (m_SelectedImage.IsNotNull()) && (m_ImageObserverTag >= 0) ) { m_SelectedImage->RemoveObserver( m_ImageObserverTag ); m_ImageObserverTag = -1; } *m_SelectedImageNodes = node; m_SelectedImage = selectedImage; m_IsTensorImage = false; m_FAImage = nullptr; m_RDImage = nullptr; m_ADImage = nullptr; m_MDImage = nullptr; m_CAImage = nullptr; m_DirectionComp1Image = nullptr; m_DirectionComp2Image = nullptr; m_AngularErrorImage = nullptr; // Add change listeners to selected objects m_ImageObserverTag = m_SelectedImage->AddObserver( itk::ModifiedEvent(), changeListener ); m_Controls->m_SimilarAnglesFrame->setVisible(false); m_Controls->m_SimilarAnglesLabel->setVisible(false); m_Controls->m_SelectedImageLabel->setText( m_SelectedImageNodes->GetNode()->GetName().c_str() ); } } //planar if(isPlanar) { if(selectedPlanar != m_SelectedPlanarFigure.GetPointer()) { MITK_DEBUG << "Planar selection changed"; somethingChanged = true; // Possibly previous change listeners if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_PlanarFigureObserverTag >= 0) ) { m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); m_PlanarFigureObserverTag = -1; } if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_InitializedObserverTag >= 0) ) { m_SelectedPlanarFigure->RemoveObserver( m_InitializedObserverTag ); m_InitializedObserverTag = -1; } m_SelectedPlanarFigure = selectedPlanar; *m_SelectedPlanarFigureNodes = node; m_CurrentFigureNodeInitialized = selectedPlanar->IsPlaced(); m_SelectedMaskNode = nullptr; m_SelectedImageMask = nullptr; m_PlanarFigureObserverTag = m_SelectedPlanarFigure->AddObserver( mitk::EndInteractionPlanarFigureEvent(), changeListener ); if(!m_CurrentFigureNodeInitialized) { typedef itk::SimpleMemberCommand< QmitkPartialVolumeAnalysisView > ITKCommandType; ITKCommandType::Pointer initializationCommand; initializationCommand = ITKCommandType::New(); // set the callback function of the member command initializationCommand->SetCallbackFunction( this, &QmitkPartialVolumeAnalysisView::PlanarFigureInitialized ); // add an observer m_InitializedObserverTag = selectedPlanar->AddObserver( mitk::EndPlacementPlanarFigureEvent(), initializationCommand ); } m_Controls->m_SelectedMaskLabel->setText( m_SelectedPlanarFigureNodes->GetNode()->GetName().c_str() ); PlanarFigureFocus(node); } } //mask this->m_Controls->m_EstimateCircle->setEnabled(isMask && selectedImage->GetDimension()==3); if(isMask && selectedImage->GetDimension()==3) { if(selectedImage != m_SelectedImage.GetPointer()) { somethingChanged = true; if ( (m_SelectedImageMask.IsNotNull()) && (m_ImageMaskObserverTag >= 0) ) { m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); m_ImageMaskObserverTag = -1; } m_SelectedMaskNode = node; m_SelectedImageMask = selectedImage; m_SelectedPlanarFigure = nullptr; m_SelectedPlanarFigureNodes->RemoveAllNodes(); m_ImageMaskObserverTag = m_SelectedImageMask->AddObserver( itk::ModifiedEvent(), changeListener ); m_Controls->m_SelectedMaskLabel->setText( m_SelectedMaskNode->GetName().c_str() ); } } //tensor image if(isTensorImage && selectedTensorImage->GetDimension()==3) { if(selectedImage != m_SelectedImage.GetPointer()) { somethingChanged = true; if ( (m_SelectedImage.IsNotNull()) && (m_ImageObserverTag >= 0) ) { m_SelectedImage->RemoveObserver( m_ImageObserverTag ); m_ImageObserverTag = -1; } *m_SelectedImageNodes = node; m_SelectedImage = selectedImage; m_IsTensorImage = true; ExtractTensorImages(selectedImage); // Add change listeners to selected objects m_ImageObserverTag = m_SelectedImage->AddObserver( itk::ModifiedEvent(), changeListener ); m_Controls->m_SimilarAnglesFrame->setVisible(true); m_Controls->m_SimilarAnglesLabel->setVisible(true); m_Controls->m_SelectedImageLabel->setText( m_SelectedImageNodes->GetNode()->GetName().c_str() ); } } } if(somethingChanged) { this->SetMeasurementInfoToRenderWindow(""); if(m_SelectedPlanarFigure.IsNull() && m_SelectedImageMask.IsNull() ) { m_Controls->m_SelectedMaskLabel->setText("mandatory"); m_Controls->m_ResampleOptionsFrame->setEnabled(false); m_Controls->m_HistogramWidget->setEnabled(false); m_Controls->m_ClassSelector->setEnabled(false); m_Controls->m_DisplayHistogramCheckbox->setEnabled(false); m_Controls->m_AdvancedCheckbox->setEnabled(false); m_Controls->frame_7->setEnabled(false); } else { m_Controls->m_ResampleOptionsFrame->setEnabled(true); m_Controls->m_HistogramWidget->setEnabled(true); m_Controls->m_ClassSelector->setEnabled(true); m_Controls->m_DisplayHistogramCheckbox->setEnabled(true); m_Controls->m_AdvancedCheckbox->setEnabled(true); m_Controls->frame_7->setEnabled(true); } // Clear statistics / histogram GUI if nothing is selected if ( m_SelectedImage.IsNull() ) { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); m_Controls->m_OpacityFrame->setEnabled(false); m_Controls->m_SelectedImageLabel->setText("mandatory"); } else { m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); m_Controls->m_OpacityFrame->setEnabled(true); } if( !m_Visible || m_SelectedImage.IsNull() || (m_SelectedPlanarFigure.IsNull() && m_SelectedImageMask.IsNull()) ) { m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_HistogramWidget->ClearItemModel(); m_CurrentStatisticsValid = false; } else { m_Controls->m_InputData->setTitle("Input Data"); this->RequestStatisticsUpdate(); } } } void QmitkPartialVolumeAnalysisView::ShowClusteringResults() { typedef itk::Image MaskImageType; mitk::Image::Pointer mask = 0; MaskImageType::Pointer itkmask = 0; if(m_IsTensorImage && m_Controls->m_SimilarAnglesSlider->value() != 0) { typedef itk::Image AngularErrorImageType; typedef mitk::ImageToItk CastType; CastType::Pointer caster = CastType::New(); caster->SetInput(m_AngularErrorImage); caster->Update(); typedef itk::BinaryThresholdImageFilter< AngularErrorImageType, MaskImageType > ThreshType; ThreshType::Pointer thresh = ThreshType::New(); thresh->SetUpperThreshold((90-m_Controls->m_SimilarAnglesSlider->value())*(itk::Math::pi/180.0)); thresh->SetInsideValue(1.0); thresh->SetInput(caster->GetOutput()); thresh->Update(); itkmask = thresh->GetOutput(); mask = mitk::Image::New(); mask->InitializeByItk(itkmask.GetPointer()); mask->SetVolume(itkmask->GetBufferPointer()); // GetDataStorage()->Remove(m_newnode); // m_newnode = mitk::DataNode::New(); // m_newnode->SetData(mask); // m_newnode->SetName("masking node"); // m_newnode->SetIntProperty( "layer", 1002 ); // GetDataStorage()->Add(m_newnode, m_SelectedImageNodes->GetNode()); } mitk::Image::Pointer clusteredImage; ClusteringType::Pointer clusterer = ClusteringType::New(); if(m_QuantifyClass==3) { if(m_IsTensorImage) { double *green_fa, *green_rd, *green_ad, *green_md; //double *greengray_fa, *greengray_rd, *greengray_ad, *greengray_md; double *gray_fa, *gray_rd, *gray_ad, *gray_md; //double *redgray_fa, *redgray_rd, *redgray_ad, *redgray_md; double *red_fa, *red_rd, *red_ad, *red_md; mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(0); mitk::Image::ConstPointer imgToCluster = tmpImg; red_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->r, mask); green_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->g, mask); gray_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->b, mask); tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(3); mitk::Image::ConstPointer imgToCluster3 = tmpImg; red_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentRGBClusteringResults->rgbChannels->r, mask); green_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentRGBClusteringResults->rgbChannels->g, mask); gray_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentRGBClusteringResults->rgbChannels->b, mask); tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(4); mitk::Image::ConstPointer imgToCluster4 = tmpImg; red_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentRGBClusteringResults->rgbChannels->r, mask); green_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentRGBClusteringResults->rgbChannels->g, mask); gray_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentRGBClusteringResults->rgbChannels->b, mask); tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(5); mitk::Image::ConstPointer imgToCluster5 = tmpImg; red_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentRGBClusteringResults->rgbChannels->r, mask); green_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentRGBClusteringResults->rgbChannels->g, mask); gray_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentRGBClusteringResults->rgbChannels->b, mask); // clipboard QString clipboardText("FA\t%1\t%2\t\t%3\t%4\t\t%5\t%6\t"); clipboardText = clipboardText .arg(red_fa[0]).arg(red_fa[1]) .arg(gray_fa[0]).arg(gray_fa[1]) .arg(green_fa[0]).arg(green_fa[1]); QString clipboardText3("RD\t%1\t%2\t\t%3\t%4\t\t%5\t%6\t"); clipboardText3 = clipboardText3 .arg(red_rd[0]).arg(red_rd[1]) .arg(gray_rd[0]).arg(gray_rd[1]) .arg(green_rd[0]).arg(green_rd[1]); QString clipboardText4("AD\t%1\t%2\t\t%3\t%4\t\t%5\t%6\t"); clipboardText4 = clipboardText4 .arg(red_ad[0]).arg(red_ad[1]) .arg(gray_ad[0]).arg(gray_ad[1]) .arg(green_ad[0]).arg(green_ad[1]); QString clipboardText5("MD\t%1\t%2\t\t%3\t%4\t\t%5\t%6"); clipboardText5 = clipboardText5 .arg(red_md[0]).arg(red_md[1]) .arg(gray_md[0]).arg(gray_md[1]) .arg(green_md[0]).arg(green_md[1]); QApplication::clipboard()->setText(clipboardText+clipboardText3+clipboardText4+clipboardText5, QClipboard::Clipboard); // now paint infos also on renderwindow QString plainInfoText("%1 %2 %3 \n"); plainInfoText = plainInfoText .arg("Red ", 20) .arg("Gray ", 20) .arg("Green", 20); QString plainInfoText0("FA:%1 ± %2%3 ± %4%5 ± %6\n"); plainInfoText0 = plainInfoText0 .arg(red_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(red_fa[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(gray_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(gray_fa[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(green_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(green_fa[1], -10, 'g', 2, QLatin1Char( ' ' )); QString plainInfoText3("RDx10³:%1 ± %2%3 ± %4%5 ± %6\n"); plainInfoText3 = plainInfoText3 .arg(1000.0 * red_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_rd[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(1000.0 * gray_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * gray_rd[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(1000.0 * green_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * green_rd[1], -10, 'g', 2, QLatin1Char( ' ' )); QString plainInfoText4("ADx10³:%1 ± %2%3 ± %4%5 ± %6\n"); plainInfoText4 = plainInfoText4 .arg(1000.0 * red_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_ad[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(1000.0 * gray_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * gray_ad[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(1000.0 * green_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * green_ad[1], -10, 'g', 2, QLatin1Char( ' ' )); QString plainInfoText5("MDx10³:%1 ± %2%3 ± %4%5 ± %6"); plainInfoText5 = plainInfoText5 .arg(1000.0 * red_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_md[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(1000.0 * gray_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * gray_md[1], -10, 'g', 2, QLatin1Char( ' ' )) .arg(1000.0 * green_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * green_md[1], -10, 'g', 2, QLatin1Char( ' ' )); this->SetMeasurementInfoToRenderWindow(plainInfoText+plainInfoText0+plainInfoText3+plainInfoText4+plainInfoText5); } else { double* green; double* gray; double* red; mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalImage(); mitk::Image::ConstPointer imgToCluster = tmpImg; red = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->r); green = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->g); gray = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->b); // clipboard QString clipboardText("%1\t%2\t\t%3\t%4\t\t%5\t%6"); clipboardText = clipboardText.arg(red[0]).arg(red[1]) .arg(gray[0]).arg(gray[1]) .arg(green[0]).arg(green[1]); QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); // now paint infos also on renderwindow QString plainInfoText("Red: %1 ± %2\nGray: %3 ± %4\nGreen: %5 ± %6"); plainInfoText = plainInfoText.arg(red[0]).arg(red[1]) .arg(gray[0]).arg(gray[1]) .arg(green[0]).arg(green[1]); this->SetMeasurementInfoToRenderWindow(plainInfoText); } clusteredImage = m_CurrentRGBClusteringResults->rgb; } else { if(m_IsTensorImage) { double *red_fa, *red_rd, *red_ad, *red_md; mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(0); mitk::Image::ConstPointer imgToCluster = tmpImg; red_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentPerformClusteringResults->clusteredImage, mask); tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(3); mitk::Image::ConstPointer imgToCluster3 = tmpImg; red_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentPerformClusteringResults->clusteredImage, mask); tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(4); mitk::Image::ConstPointer imgToCluster4 = tmpImg; red_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentPerformClusteringResults->clusteredImage, mask); tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(5); mitk::Image::ConstPointer imgToCluster5 = tmpImg; red_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentPerformClusteringResults->clusteredImage, mask); // clipboard QString clipboardText("FA\t%1\t%2\t"); clipboardText = clipboardText .arg(red_fa[0]).arg(red_fa[1]); QString clipboardText3("RD\t%1\t%2\t"); clipboardText3 = clipboardText3 .arg(red_rd[0]).arg(red_rd[1]); QString clipboardText4("AD\t%1\t%2\t"); clipboardText4 = clipboardText4 .arg(red_ad[0]).arg(red_ad[1]); QString clipboardText5("MD\t%1\t%2\t"); clipboardText5 = clipboardText5 .arg(red_md[0]).arg(red_md[1]); QApplication::clipboard()->setText(clipboardText+clipboardText3+clipboardText4+clipboardText5, QClipboard::Clipboard); // now paint infos also on renderwindow QString plainInfoText("%1 \n"); plainInfoText = plainInfoText .arg("Red ", 20); QString plainInfoText0("FA:%1 ± %2\n"); plainInfoText0 = plainInfoText0 .arg(red_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(red_fa[1], -10, 'g', 2, QLatin1Char( ' ' )); QString plainInfoText3("RDx10³:%1 ± %2\n"); plainInfoText3 = plainInfoText3 .arg(1000.0 * red_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_rd[1], -10, 'g', 2, QLatin1Char( ' ' )); QString plainInfoText4("ADx10³:%1 ± %2\n"); plainInfoText4 = plainInfoText4 .arg(1000.0 * red_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_ad[1], -10, 'g', 2, QLatin1Char( ' ' )); QString plainInfoText5("MDx10³:%1 ± %2"); plainInfoText5 = plainInfoText5 .arg(1000.0 * red_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_md[1], -10, 'g', 2, QLatin1Char( ' ' )); this->SetMeasurementInfoToRenderWindow(plainInfoText+plainInfoText0+plainInfoText3+plainInfoText4+plainInfoText5); } else { double* quant; mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalImage(); mitk::Image::ConstPointer imgToCluster = tmpImg; quant = clusterer->PerformQuantification(imgToCluster, m_CurrentPerformClusteringResults->clusteredImage); // clipboard QString clipboardText("%1\t%2"); clipboardText = clipboardText.arg(quant[0]).arg(quant[1]); QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); // now paint infos also on renderwindow QString plainInfoText("Measurement: %1 ± %2"); plainInfoText = plainInfoText.arg(quant[0]).arg(quant[1]); this->SetMeasurementInfoToRenderWindow(plainInfoText); } clusteredImage = m_CurrentPerformClusteringResults->displayImage; } if(mask.IsNotNull()) { typedef itk::Image,3> RGBImageType; typedef mitk::ImageToItk ClusterCasterType; ClusterCasterType::Pointer clCaster = ClusterCasterType::New(); clCaster->SetInput(clusteredImage); clCaster->Update(); clCaster->GetOutput(); typedef itk::MaskImageFilter< RGBImageType, MaskImageType, RGBImageType > MaskType; MaskType::Pointer masker = MaskType::New(); masker->SetInput1(clCaster->GetOutput()); masker->SetInput2(itkmask); masker->Update(); clusteredImage = mitk::Image::New(); clusteredImage->InitializeByItk(masker->GetOutput()); clusteredImage->SetVolume(masker->GetOutput()->GetBufferPointer()); } if(m_ClusteringResult.IsNotNull()) { this->GetDataStorage()->Remove(m_ClusteringResult); } m_ClusteringResult = mitk::DataNode::New(); m_ClusteringResult->SetBoolProperty("helper object", true); m_ClusteringResult->SetIntProperty( "layer", 1000 ); m_ClusteringResult->SetBoolProperty("texture interpolation", m_TexIsOn); m_ClusteringResult->SetData(clusteredImage); m_ClusteringResult->SetName("Clusterprobs"); this->GetDataStorage()->Add(m_ClusteringResult, m_SelectedImageNodes->GetNode()); if(m_SelectedPlanarFigure.IsNotNull() && m_SelectedPlanarFigureNodes->GetNode().IsNotNull()) { m_SelectedPlanarFigureNodes->GetNode()->SetIntProperty( "layer", 1001 ); } this->RequestRenderWindowUpdate(); } void QmitkPartialVolumeAnalysisView::UpdateStatistics() { if(!m_CurrentFigureNodeInitialized && m_SelectedPlanarFigure.IsNotNull()) { MITK_DEBUG << "Selected planar figure not initialized. No stats calculation performed."; return; } // Remove any cached images that are no longer referenced elsewhere this->RemoveOrphanImages(); if ( m_SelectedImage.IsNotNull() ) { // Check if a the selected image is a multi-channel image. If yes, statistics // cannot be calculated currently. if ( !m_IsTensorImage && m_SelectedImage->GetPixelType().GetNumberOfComponents() > 1 ) { QMessageBox::information( nullptr, "Warning", "Non-tensor multi-component images not supported."); m_Controls->m_HistogramWidget->ClearItemModel(); m_CurrentStatisticsValid = false; return; } m_CurrentStatisticsCalculator = nullptr; if(!m_IsTensorImage) { // Retrieve HistogramStatisticsCalculator from has map (or create a new one // for this image if non-existant) PartialVolumeAnalysisMapType::iterator it = m_PartialVolumeAnalysisMap.find( m_SelectedImage ); if ( it != m_PartialVolumeAnalysisMap.end() ) { m_CurrentStatisticsCalculator = it->second; } } if(m_CurrentStatisticsCalculator.IsNull()) { m_CurrentStatisticsCalculator = mitk::PartialVolumeAnalysisHistogramCalculator::New(); m_CurrentStatisticsCalculator->SetPlanarFigureThickness(m_Controls->m_PlanarFiguresThickness->value()); if(m_IsTensorImage) { m_CurrentStatisticsCalculator->SetImage( m_CAImage ); m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_FAImage ); m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_DirectionComp1Image ); m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_DirectionComp2Image ); m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_RDImage ); m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_ADImage ); m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_MDImage ); } else { m_CurrentStatisticsCalculator->SetImage( m_SelectedImage ); } m_PartialVolumeAnalysisMap[m_SelectedImage] = m_CurrentStatisticsCalculator; MITK_DEBUG << "Creating StatisticsCalculator"; } std::string maskName; std::string maskType; unsigned int maskDimension; if ( m_SelectedImageMask.IsNotNull() ) { mitk::PixelType pixelType = m_SelectedImageMask->GetPixelType(); MITK_DEBUG << pixelType.GetPixelTypeAsString(); if(pixelType.GetComponentTypeAsString() == "char") { MITK_DEBUG << "Pixel type is char instead of uchar"; return; } if(pixelType.GetBitsPerComponent() == 16) { //convert from ushort to uchar typedef itk::Image UCharImageType; UCharImageType::Pointer charImage; if(pixelType.GetComponentTypeAsString() == "short" ) { typedef itk::Image ShortImageType; ShortImageType::Pointer shortImage; mitk::CastToItkImage( m_SelectedImageMask, shortImage ); typedef itk::ImageDuplicator< ShortImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( shortImage ); duplicator->Update(); typedef itk::CastImageFilter ImageCasterType; ImageCasterType::Pointer caster = ImageCasterType::New(); caster->SetInput( duplicator->GetOutput() ); caster->Update(); charImage = caster->GetOutput(); } else { typedef itk::Image UShortImageType; UShortImageType::Pointer shortImage; mitk::CastToItkImage( m_SelectedImageMask, shortImage ); typedef itk::ImageDuplicator< UShortImageType > DuplicatorType; DuplicatorType::Pointer duplicator = DuplicatorType::New(); duplicator->SetInputImage( shortImage ); duplicator->Update(); typedef itk::CastImageFilter ImageCasterType; ImageCasterType::Pointer caster = ImageCasterType::New(); caster->SetInput( duplicator->GetOutput() ); caster->Update(); charImage = caster->GetOutput(); } m_SelectedImageMask = nullptr; m_SelectedImageMask = mitk::Image::New(); m_SelectedImageMask->InitializeByItk( charImage.GetPointer() ); m_SelectedImageMask->SetVolume( charImage->GetBufferPointer() ); mitk::CastToMitkImage(charImage, m_SelectedImageMask); } m_CurrentStatisticsCalculator->SetImageMask( m_SelectedImageMask ); m_CurrentStatisticsCalculator->SetMaskingModeToImage(); maskName = m_SelectedMaskNode->GetName(); maskType = m_SelectedImageMask->GetNameOfClass(); maskDimension = 3; std::stringstream maskLabel; maskLabel << maskName; if ( maskDimension > 0 ) { maskLabel << " [" << maskDimension << "D " << maskType << "]"; } m_Controls->m_SelectedMaskLabel->setText( maskLabel.str().c_str() ); } else if ( m_SelectedPlanarFigure.IsNotNull() && m_SelectedPlanarFigureNodes->GetNode().IsNotNull()) { m_CurrentStatisticsCalculator->SetPlanarFigure( m_SelectedPlanarFigure ); m_CurrentStatisticsCalculator->SetMaskingModeToPlanarFigure(); maskName = m_SelectedPlanarFigureNodes->GetNode()->GetName(); maskType = m_SelectedPlanarFigure->GetNameOfClass(); maskDimension = 2; } else { m_CurrentStatisticsCalculator->SetMaskingModeToNone(); maskName = "-"; maskType = ""; maskDimension = 0; } bool statisticsChanged = false; bool statisticsCalculationSuccessful = false; // Initialize progress bar mitk::ProgressBar::GetInstance()->AddStepsToDo( 100 ); // Install listener for progress events and initialize progress bar typedef itk::SimpleMemberCommand< QmitkPartialVolumeAnalysisView > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &QmitkPartialVolumeAnalysisView::UpdateProgressBar ); unsigned long progressObserverTag = m_CurrentStatisticsCalculator ->AddObserver( itk::ProgressEvent(), progressListener ); ClusteringType::ParamsType *cparams = 0; ClusteringType::ClusterResultType *cresult = 0; ClusteringType::HistType *chist = 0; try { m_CurrentStatisticsCalculator->SetNumberOfBins(m_Controls->m_NumberBins->text().toInt()); m_CurrentStatisticsCalculator->SetUpsamplingFactor(m_Controls->m_Upsampling->text().toDouble()); m_CurrentStatisticsCalculator->SetGaussianSigma(m_Controls->m_GaussianSigma->text().toDouble()); // Compute statistics statisticsChanged = m_CurrentStatisticsCalculator->ComputeStatistics( ); mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalImage(); mitk::Image::ConstPointer imgToCluster = tmpImg; if(imgToCluster.IsNotNull()) { // perform clustering const HistogramType *histogram = m_CurrentStatisticsCalculator->GetHistogram( ); if(histogram != nullptr) { ClusteringType::Pointer clusterer = ClusteringType::New(); clusterer->SetStepsNumIntegration(200); clusterer->SetMaxIt(1000); mitk::Image::Pointer pFiberImg; if(m_QuantifyClass==3) { if(m_Controls->m_Quantiles->isChecked()) { m_CurrentRGBClusteringResults = clusterer->PerformRGBQuantiles(imgToCluster, histogram, m_Controls->m_q1->value(),m_Controls->m_q2->value()); } else { m_CurrentRGBClusteringResults = clusterer->PerformRGBClustering(imgToCluster, histogram); } pFiberImg = m_CurrentRGBClusteringResults->rgbChannels->r; cparams = m_CurrentRGBClusteringResults->params; cresult = m_CurrentRGBClusteringResults->result; chist = m_CurrentRGBClusteringResults->hist; } else { if(m_Controls->m_Quantiles->isChecked()) { m_CurrentPerformClusteringResults = clusterer->PerformQuantiles(imgToCluster, histogram, m_Controls->m_q1->value(),m_Controls->m_q2->value()); } else { m_CurrentPerformClusteringResults = clusterer->PerformClustering(imgToCluster, histogram, m_QuantifyClass); } pFiberImg = m_CurrentPerformClusteringResults->clusteredImage; cparams = m_CurrentPerformClusteringResults->params; cresult = m_CurrentPerformClusteringResults->result; chist = m_CurrentPerformClusteringResults->hist; } if(m_IsTensorImage) { m_AngularErrorImage = clusterer->CaculateAngularErrorImage( m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(1), m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(2), pFiberImg); // GetDataStorage()->Remove(m_newnode2); // m_newnode2 = mitk::DataNode::New(); // m_newnode2->SetData(m_AngularErrorImage); // m_newnode2->SetName(("AngularError")); // m_newnode2->SetIntProperty( "layer", 1003 ); // GetDataStorage()->Add(m_newnode2, m_SelectedImageNodes->GetNode()); // newnode = mitk::DataNode::New(); // newnode->SetData(m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(1)); // newnode->SetName(("Comp1")); // GetDataStorage()->Add(newnode, m_SelectedImageNodes->GetNode()); // newnode = mitk::DataNode::New(); // newnode->SetData(m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(2)); // newnode->SetName(("Comp2")); // GetDataStorage()->Add(newnode, m_SelectedImageNodes->GetNode()); } ShowClusteringResults(); } } statisticsCalculationSuccessful = true; } catch ( const std::runtime_error &e ) { QMessageBox::information( nullptr, "Warning", e.what()); } catch ( const std::exception &e ) { MITK_ERROR << "Caught exception: " << e.what(); QMessageBox::information( nullptr, "Warning", e.what()); } m_CurrentStatisticsCalculator->RemoveObserver( progressObserverTag ); // Make sure that progress bar closes mitk::ProgressBar::GetInstance()->Progress( 100 ); if ( statisticsCalculationSuccessful ) { if ( statisticsChanged ) { // Do not show any error messages m_CurrentStatisticsValid = true; } // m_Controls->m_HistogramWidget->SetHistogramModeToDirectHistogram(); m_Controls->m_HistogramWidget->SetParameters( cparams, cresult, chist ); // m_Controls->m_HistogramWidget->UpdateItemModelFromHistogram(); } else { m_Controls->m_SelectedMaskLabel->setText("mandatory"); // Clear statistics and histogram m_Controls->m_HistogramWidget->ClearItemModel(); m_CurrentStatisticsValid = false; // If a (non-closed) PlanarFigure is selected, display a line profile widget if ( m_SelectedPlanarFigure.IsNotNull() ) { // TODO: enable line profile widget //m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 1 ); //m_Controls->m_LineProfileWidget->SetImage( m_SelectedImage ); //m_Controls->m_LineProfileWidget->SetPlanarFigure( m_SelectedPlanarFigure ); //m_Controls->m_LineProfileWidget->UpdateItemModelFromPath(); } } } } void QmitkPartialVolumeAnalysisView::SetMeasurementInfoToRenderWindow(const QString& text) { FindRenderWindow(m_SelectedPlanarFigureNodes->GetNode()); if(m_LastRenderWindow != m_SelectedRenderWindow) { if(m_LastRenderWindow) { QObject::disconnect( m_LastRenderWindow, SIGNAL( destroyed(QObject*) ) , this, SLOT( OnRenderWindowDelete(QObject*) ) ); } m_LastRenderWindow = m_SelectedRenderWindow; if(m_LastRenderWindow) { QObject::connect( m_LastRenderWindow, SIGNAL( destroyed(QObject*) ) , this, SLOT( OnRenderWindowDelete(QObject*) ) ); } } if(m_LastRenderWindow && m_SelectedPlanarFigureNodes->GetNode().IsNotNull()) { if (!text.isEmpty()) { m_MeasurementInfoAnnotation->SetText(1, text.toLatin1().data()); - mitk::VtkLayerController::GetInstance(m_LastRenderWindow->GetRenderWindow())->InsertForegroundRenderer( + mitk::VtkLayerController::GetInstance(m_LastRenderWindow->GetVtkRenderWindow())->InsertForegroundRenderer( m_MeasurementInfoRenderer, true); } else { if (mitk::VtkLayerController::GetInstance( - m_LastRenderWindow->GetRenderWindow()) ->IsRendererInserted( + m_LastRenderWindow->GetVtkRenderWindow()) ->IsRendererInserted( m_MeasurementInfoRenderer)) - mitk::VtkLayerController::GetInstance(m_LastRenderWindow->GetRenderWindow())->RemoveRenderer( + mitk::VtkLayerController::GetInstance(m_LastRenderWindow->GetVtkRenderWindow())->RemoveRenderer( m_MeasurementInfoRenderer); } } else { mitk::IRenderWindowPart* renderWindowPart = this->GetRenderWindowPart(); if ( renderWindowPart == nullptr ) { return; } if (!text.isEmpty()) { m_MeasurementInfoAnnotation->SetText(1, text.toLatin1().data()); - mitk::VtkLayerController::GetInstance(renderWindowPart->GetQmitkRenderWindow("axial")->GetRenderWindow())->InsertForegroundRenderer( + mitk::VtkLayerController::GetInstance(renderWindowPart->GetQmitkRenderWindow("axial")->GetVtkRenderWindow())->InsertForegroundRenderer( m_MeasurementInfoRenderer, true); } else { if (mitk::VtkLayerController::GetInstance( - renderWindowPart->GetQmitkRenderWindow("axial")->GetRenderWindow()) ->IsRendererInserted( + renderWindowPart->GetQmitkRenderWindow("axial")->GetVtkRenderWindow()) ->IsRendererInserted( m_MeasurementInfoRenderer)) - mitk::VtkLayerController::GetInstance(renderWindowPart->GetQmitkRenderWindow("axial")->GetRenderWindow())->RemoveRenderer( + mitk::VtkLayerController::GetInstance(renderWindowPart->GetQmitkRenderWindow("axial")->GetVtkRenderWindow())->RemoveRenderer( m_MeasurementInfoRenderer); } } } void QmitkPartialVolumeAnalysisView::UpdateProgressBar() { mitk::ProgressBar::GetInstance()->Progress(); } void QmitkPartialVolumeAnalysisView::RequestStatisticsUpdate() { if ( !m_StatisticsUpdatePending ) { QApplication::postEvent( this, new QmitkRequestStatisticsUpdateEvent ); m_StatisticsUpdatePending = true; } } void QmitkPartialVolumeAnalysisView::RemoveOrphanImages() { PartialVolumeAnalysisMapType::iterator it = m_PartialVolumeAnalysisMap.begin(); while ( it != m_PartialVolumeAnalysisMap.end() ) { mitk::Image *image = it->first; mitk::PartialVolumeAnalysisHistogramCalculator *calculator = it->second; ++it; mitk::NodePredicateData::Pointer hasImage = mitk::NodePredicateData::New( image ); if ( this->GetDataStorage()->GetNode( hasImage ) == nullptr ) { if ( m_SelectedImage == image ) { m_SelectedImage = nullptr; m_SelectedImageNodes->RemoveAllNodes(); } if ( m_CurrentStatisticsCalculator == calculator ) { m_CurrentStatisticsCalculator = nullptr; } m_PartialVolumeAnalysisMap.erase( image ); it = m_PartialVolumeAnalysisMap.begin(); } } } void QmitkPartialVolumeAnalysisView::ExtractTensorImages( mitk::Image::Pointer tensorimage) { typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; typedef mitk::ImageToItk CastType; CastType::Pointer caster = CastType::New(); caster->SetInput(tensorimage); caster->Update(); TensorImageType::Pointer image = caster->GetOutput(); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(image ); measurementsCalculator->SetMeasure(MeasurementsType::FA); measurementsCalculator->Update(); MeasurementsType::OutputImageType::Pointer fa = measurementsCalculator->GetOutput(); m_FAImage = mitk::Image::New(); m_FAImage->InitializeByItk(fa.GetPointer()); m_FAImage->SetVolume(fa->GetBufferPointer()); // mitk::DataNode::Pointer node = mitk::DataNode::New(); // node->SetData(m_FAImage); // GetDataStorage()->Add(node); measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(image ); measurementsCalculator->SetMeasure(MeasurementsType::CA); measurementsCalculator->Update(); MeasurementsType::OutputImageType::Pointer ca = measurementsCalculator->GetOutput(); m_CAImage = mitk::Image::New(); m_CAImage->InitializeByItk(ca.GetPointer()); m_CAImage->SetVolume(ca->GetBufferPointer()); // node = mitk::DataNode::New(); // node->SetData(m_CAImage); // GetDataStorage()->Add(node); measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(image ); measurementsCalculator->SetMeasure(MeasurementsType::RD); measurementsCalculator->Update(); MeasurementsType::OutputImageType::Pointer rd = measurementsCalculator->GetOutput(); m_RDImage = mitk::Image::New(); m_RDImage->InitializeByItk(rd.GetPointer()); m_RDImage->SetVolume(rd->GetBufferPointer()); // node = mitk::DataNode::New(); // node->SetData(m_CAImage); // GetDataStorage()->Add(node); measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(image ); measurementsCalculator->SetMeasure(MeasurementsType::AD); measurementsCalculator->Update(); MeasurementsType::OutputImageType::Pointer ad = measurementsCalculator->GetOutput(); m_ADImage = mitk::Image::New(); m_ADImage->InitializeByItk(ad.GetPointer()); m_ADImage->SetVolume(ad->GetBufferPointer()); // node = mitk::DataNode::New(); // node->SetData(m_CAImage); // GetDataStorage()->Add(node); measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(image ); measurementsCalculator->SetMeasure(MeasurementsType::RA); measurementsCalculator->Update(); MeasurementsType::OutputImageType::Pointer md = measurementsCalculator->GetOutput(); m_MDImage = mitk::Image::New(); m_MDImage->InitializeByItk(md.GetPointer()); m_MDImage->SetVolume(md->GetBufferPointer()); // node = mitk::DataNode::New(); // node->SetData(m_CAImage); // GetDataStorage()->Add(node); typedef DirectionsFilterType::PeakImageType DirImageType; DirectionsFilterType::Pointer dirFilter = DirectionsFilterType::New(); dirFilter->SetInput(image ); dirFilter->SetUsePolarCoordinates(true); dirFilter->Update(); ItkUcharImgType::Pointer numDirImage = dirFilter->GetOutput(); DirImageType::Pointer dirImage = dirFilter->GetPeakImage(); itk::ImageRegionIterator itd(dirImage, dirImage->GetLargestPossibleRegion()); itd.GoToBegin(); while( !itd.IsAtEnd() ) { DirImageType::PixelType direction = itd.Get(); direction = fabs(direction); itd.Set(direction); ++itd; } typedef itk::Image CompImageType; CompImageType::Pointer comp1 = CompImageType::New(); comp1->SetSpacing( numDirImage->GetSpacing() ); // Set the image spacing comp1->SetOrigin( numDirImage->GetOrigin() ); // Set the image origin comp1->SetDirection( numDirImage->GetDirection() ); // Set the image direction comp1->SetRegions( numDirImage->GetLargestPossibleRegion() ); comp1->Allocate(); CompImageType::Pointer comp2 = CompImageType::New(); comp2->SetSpacing( numDirImage->GetSpacing() ); // Set the image spacing comp2->SetOrigin( numDirImage->GetOrigin() ); // Set the image origin comp2->SetDirection( numDirImage->GetDirection() ); // Set the image direction comp2->SetRegions( numDirImage->GetLargestPossibleRegion() ); comp2->Allocate(); itk::ImageRegionIterator it1(comp1, comp1->GetLargestPossibleRegion()); itk::ImageRegionIterator it2(comp2, comp2->GetLargestPossibleRegion()); it1.GoToBegin(); it2.GoToBegin(); while( !it2.IsAtEnd() ) { DirImageType::IndexType peakIndex; peakIndex[0] = it1.GetIndex()[0]; peakIndex[1] = it1.GetIndex()[1]; peakIndex[2] = it1.GetIndex()[2]; peakIndex[3] = 1; it1.Set(dirImage->GetPixel(peakIndex)); peakIndex[3] = 2; it2.Set(dirImage->GetPixel(peakIndex)); ++it1; ++it2; } m_DirectionComp1Image = mitk::Image::New(); m_DirectionComp1Image->InitializeByItk(comp1.GetPointer()); m_DirectionComp1Image->SetVolume(comp1->GetBufferPointer()); m_DirectionComp2Image = mitk::Image::New(); m_DirectionComp2Image->InitializeByItk(comp2.GetPointer()); m_DirectionComp2Image->SetVolume(comp2->GetBufferPointer()); } void QmitkPartialVolumeAnalysisView::OnRenderWindowDelete(QObject * obj) { if(obj == m_LastRenderWindow) m_LastRenderWindow = 0; if(obj == m_SelectedRenderWindow) m_SelectedRenderWindow = 0; } bool QmitkPartialVolumeAnalysisView::event( QEvent *event ) { if ( event->type() == (QEvent::Type) QmitkRequestStatisticsUpdateEvent::StatisticsUpdateRequest ) { // Update statistics m_StatisticsUpdatePending = false; this->UpdateStatistics(); return true; } return false; } void QmitkPartialVolumeAnalysisView::Activated() { mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDataStorage()->GetAll(); mitk::DataNode* node = 0; mitk::PlanarFigure* figure = 0; mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; // finally add all nodes to the model for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() ; it++) { node = const_cast(it->Value().GetPointer()); figure = dynamic_cast(node->GetData()); if(figure) { figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( node ); } } } } void QmitkPartialVolumeAnalysisView::Deactivated() { } void QmitkPartialVolumeAnalysisView::ActivatedZombieView(berry::IWorkbenchPartReference::Pointer) { this->SetMeasurementInfoToRenderWindow(""); mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDataStorage()->GetAll(); mitk::DataNode* node = 0; mitk::PlanarFigure* figure = 0; mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; // finally add all nodes to the model for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() ; it++) { node = const_cast(it->Value().GetPointer()); figure = dynamic_cast(node->GetData()); if(figure) { figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); if(figureInteractor) figureInteractor->SetDataNode( nullptr ); } } } void QmitkPartialVolumeAnalysisView::Hidden() { if (m_ClusteringResult.IsNotNull()) { this->GetDataStorage()->Remove(m_ClusteringResult); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } Select(nullptr, true, true); m_Visible = false; } void QmitkPartialVolumeAnalysisView::Visible() { m_Visible = true; berry::IWorkbenchPart::Pointer bla; if (!this->GetCurrentSelection().empty()) { this->OnSelectionChanged(bla, this->GetCurrentSelection()); } else { this->OnSelectionChanged(bla, this->GetDataManagerSelection()); } } void QmitkPartialVolumeAnalysisView::GreenRadio(bool checked) { if(checked) { m_Controls->m_PartialVolumeRadio->setChecked(false); m_Controls->m_BlueRadio->setChecked(false); m_Controls->m_AllRadio->setChecked(false); m_Controls->m_ExportClusteringResultsButton->setEnabled(true); } m_QuantifyClass = 0; RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::PartialVolumeRadio(bool checked) { if(checked) { m_Controls->m_GreenRadio->setChecked(false); m_Controls->m_BlueRadio->setChecked(false); m_Controls->m_AllRadio->setChecked(false); m_Controls->m_ExportClusteringResultsButton->setEnabled(true); } m_QuantifyClass = 1; RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::BlueRadio(bool checked) { if(checked) { m_Controls->m_PartialVolumeRadio->setChecked(false); m_Controls->m_GreenRadio->setChecked(false); m_Controls->m_AllRadio->setChecked(false); m_Controls->m_ExportClusteringResultsButton->setEnabled(true); } m_QuantifyClass = 2; RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::AllRadio(bool checked) { if(checked) { m_Controls->m_BlueRadio->setChecked(false); m_Controls->m_PartialVolumeRadio->setChecked(false); m_Controls->m_GreenRadio->setChecked(false); m_Controls->m_ExportClusteringResultsButton->setEnabled(false); } m_QuantifyClass = 3; RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::NumberBinsChangedSlider(int) { m_Controls->m_NumberBins->setText(QString("%1").arg(m_Controls->m_NumberBinsSlider->value()*5.0)); } void QmitkPartialVolumeAnalysisView::UpsamplingChangedSlider(int) { m_Controls->m_Upsampling->setText(QString("%1").arg(m_Controls->m_UpsamplingSlider->value()/10.0)); } void QmitkPartialVolumeAnalysisView::GaussianSigmaChangedSlider(int) { m_Controls->m_GaussianSigma->setText(QString("%1").arg(m_Controls->m_GaussianSigmaSlider->value()/100.0)); } void QmitkPartialVolumeAnalysisView::SimilarAnglesChangedSlider(int) { m_Controls->m_SimilarAngles->setText(QString("%1°").arg(90-m_Controls->m_SimilarAnglesSlider->value())); ShowClusteringResults(); } void QmitkPartialVolumeAnalysisView::OpacityChangedSlider(int v ) { if(m_SelectedImageNodes->GetNode().IsNotNull()) { float opacImag = 1.0f-(v-5)/5.0f; opacImag = opacImag < 0 ? 0 : opacImag; m_SelectedImageNodes->GetNode()->SetFloatProperty("opacity", opacImag); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } if(m_ClusteringResult.IsNotNull()) { float opacClust = v/5.0f; opacClust = opacClust > 1 ? 1 : opacClust; m_ClusteringResult->SetFloatProperty("opacity", opacClust); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } void QmitkPartialVolumeAnalysisView::NumberBinsReleasedSlider( ) { RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::UpsamplingReleasedSlider( ) { RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::GaussianSigmaReleasedSlider( ) { RequestStatisticsUpdate(); } void QmitkPartialVolumeAnalysisView::SimilarAnglesReleasedSlider( ) { } void QmitkPartialVolumeAnalysisView::ToClipBoard() { std::vector* > vals = m_Controls->m_HistogramWidget->m_Vals; QString clipboardText; for (std::vector* >::iterator it = vals.begin(); it != vals.end(); ++it) { for (std::vector::iterator it2 = (**it).begin(); it2 != (**it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); } clipboardText.append(QString("\n")); } QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); } void QmitkPartialVolumeAnalysisView::PropertyChanged(const mitk::DataNode* /*node*/, const mitk::BaseProperty* /*prop*/) { } void QmitkPartialVolumeAnalysisView::NodeChanged(const mitk::DataNode* /*node*/) { } void QmitkPartialVolumeAnalysisView::NodeRemoved(const mitk::DataNode* node) { if (dynamic_cast(node->GetData())) this->GetDataStorage()->Remove(m_ClusteringResult); if( node == m_SelectedPlanarFigureNodes->GetNode().GetPointer() || node == m_SelectedMaskNode.GetPointer() ) { this->Select(nullptr,true,false); SetMeasurementInfoToRenderWindow(""); } if( node == m_SelectedImageNodes->GetNode().GetPointer() ) { this->Select(nullptr,false,true); SetMeasurementInfoToRenderWindow(""); } } void QmitkPartialVolumeAnalysisView::NodeAddedInDataStorage(const mitk::DataNode* node) { if(!m_Visible) return; mitk::DataNode* nonConstNode = const_cast(node); mitk::PlanarFigure* figure = dynamic_cast(nonConstNode->GetData()); if(figure) { // set interactor for new node (if not already set) mitk::PlanarFigureInteractor::Pointer figureInteractor = dynamic_cast(node->GetDataInteractor().GetPointer()); if(figureInteractor.IsNull()) { figureInteractor = mitk::PlanarFigureInteractor::New(); us::Module* planarFigureModule = us::ModuleRegistry::GetModule( "MitkPlanarFigure" ); figureInteractor->LoadStateMachine("PlanarFigureInteraction.xml", planarFigureModule ); figureInteractor->SetEventConfig( "PlanarFigureConfig.xml", planarFigureModule ); figureInteractor->SetDataNode( nonConstNode ); } // remove uninitialized old planars if( m_SelectedPlanarFigureNodes->GetNode().IsNotNull() && m_CurrentFigureNodeInitialized == false ) { this->GetDataStorage()->Remove(m_SelectedPlanarFigureNodes->GetNode()); } } } void QmitkPartialVolumeAnalysisView::TextIntON() { if(m_ClusteringResult.IsNotNull()) { if(m_TexIsOn) { m_Controls->m_TextureIntON->setIcon(*m_IconTexOFF); } else { m_Controls->m_TextureIntON->setIcon(*m_IconTexON); } m_ClusteringResult->SetBoolProperty("texture interpolation", !m_TexIsOn); m_TexIsOn = !m_TexIsOn; this->RequestRenderWindowUpdate(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp index 9fe1973..149dc9e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp @@ -1,931 +1,931 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkTensorReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include #include #include #include #include // itk includes #include "itkTimeProbe.h" //#include "itkTensor.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "mitkTeemDiffusionTensor3DReconstructionImageFilter.h" #include "itkDiffusionTensor3DReconstructionImageFilter.h" #include "itkTensorImageToDiffusionImageFilter.h" #include "itkPointShell.h" #include "itkVector.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkTensorReconstructionWithEigenvalueCorrectionFilter.h" #include "mitkImageCast.h" #include "mitkImageAccessByItk.h" #include #include #include "mitkProperties.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "mitkLookupTableProperty.h" #include "mitkLookupTable.h" #include "mitkImageStatisticsHolder.h" #include #include #include #include #include #include #include #include const std::string QmitkTensorReconstructionView::VIEW_ID = "org.mitk.views.tensorreconstruction"; typedef mitk::TensorImage::ScalarPixelType TTensorPixelType; typedef mitk::TensorImage::PixelType TensorPixelType; typedef mitk::TensorImage::ItkTensorImageType TensorImageType; using namespace berry; QmitkTensorReconstructionView::QmitkTensorReconstructionView() : QmitkAbstractView(), m_Controls(nullptr) { } QmitkTensorReconstructionView::~QmitkTensorReconstructionView() { } void QmitkTensorReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkTensorReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); Advanced1CheckboxClicked(); } } void QmitkTensorReconstructionView::SetFocus() { m_Controls->m_Advanced1->setFocus(); } void QmitkTensorReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_StartReconstruction), SIGNAL(clicked()), this, SLOT(Reconstruct()) ); connect( (QObject*)(m_Controls->m_Advanced1), SIGNAL(clicked()), this, SLOT(Advanced1CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_TensorsToDWIButton), SIGNAL(clicked()), this, SLOT(TensorsToDWI()) ); connect( (QObject*)(m_Controls->m_TensorsToOdfButton), SIGNAL(clicked()), this, SLOT(TensorsToOdf()) ); connect( (QObject*)(m_Controls->m_ResidualButton), SIGNAL(clicked()), this, SLOT(ResidualCalculation()) ); connect( (QObject*)(m_Controls->m_PerSliceView), SIGNAL(pointSelected(int, int)), this, SLOT(ResidualClicked(int, int)) ); connect( (QObject*)(m_Controls->m_TensorReconstructionThreshold), SIGNAL(valueChanged(int)), this, SLOT(PreviewThreshold(int)) ); m_Controls->m_ResidualTab->setVisible(false); m_Controls->m_PercentagesOfOutliers->setVisible(false); m_Controls->m_DwiBox->SetDataStorage(this->GetDataStorage()); mitk::NodePredicateIsDWI::Pointer isDwi = mitk::NodePredicateIsDWI::New(); m_Controls->m_DwiBox->SetPredicate( isDwi ); connect( (QObject*)(m_Controls->m_DwiBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); m_Controls->m_OdfBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isDti = mitk::TNodePredicateDataType::New(); m_Controls->m_OdfBox->SetPredicate( isDti ); connect( (QObject*)(m_Controls->m_OdfBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); m_Controls->m_DtiBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_DtiBox->SetPredicate( isDti ); connect( (QObject*)(m_Controls->m_DtiBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); m_Controls->m_ResBox1->SetDataStorage(this->GetDataStorage()); m_Controls->m_ResBox1->SetPredicate( isDwi ); connect( (QObject*)(m_Controls->m_ResBox1), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); m_Controls->m_ResBox2->SetDataStorage(this->GetDataStorage()); m_Controls->m_ResBox2->SetPredicate( isDti ); connect( (QObject*)(m_Controls->m_ResBox2), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); } } void QmitkTensorReconstructionView::ResidualClicked(int slice, int volume) { mitk::Image* diffImage; mitk::DataNode::Pointer correctNode; mitk::BaseGeometry* geometry; if (m_Controls->m_ResBox1->GetSelectedNode().IsNotNull()) { diffImage = static_cast(m_Controls->m_ResBox1->GetSelectedNode()->GetData()); geometry = m_Controls->m_ResBox1->GetSelectedNode()->GetData()->GetGeometry(); // Remember the node whose display index must be updated correctNode = mitk::DataNode::New(); correctNode = m_Controls->m_ResBox1->GetSelectedNode(); GradientDirectionContainerType::ConstPointer dirs = mitk::DiffusionPropertyHelper::GetGradientContainer(diffImage); for(itk::SizeValueType i=0; iSize() && i<=(itk::SizeValueType)volume; i++) { GradientDirectionType grad = dirs->ElementAt(i); // check if image is b0 weighted if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001) { volume++; } } QString pos = "Volume: "; pos.append(QString::number(volume)); pos.append(", Slice: "); pos.append(QString::number(slice)); m_Controls->m_PositionLabel->setText(pos); if(correctNode) { int oldDisplayVal; correctNode->GetIntProperty("DisplayChannel", oldDisplayVal); QString oldVal = QString::number(oldDisplayVal); QString newVal = QString::number(volume); correctNode->SetIntProperty("DisplayChannel",volume); correctNode->SetSelected(true); this->FirePropertyChanged("DisplayChannel", oldVal, newVal); correctNode->UpdateOutputInformation(); mitk::Point3D p3 = this->GetRenderWindowPart()->GetSelectedPosition(); itk::Index<3> ix; geometry->WorldToIndex(p3, ix); // ix[2] = slice; mitk::Vector3D vec; vec[0] = ix[0]; vec[1] = ix[1]; vec[2] = slice; mitk::Vector3D v3New; geometry->IndexToWorld(vec, v3New); mitk::Point3D origin = geometry->GetOrigin(); mitk::Point3D p3New; p3New[0] = v3New[0] + origin[0]; p3New[1] = v3New[1] + origin[1]; p3New[2] = v3New[2] + origin[2]; this->GetRenderWindowPart()->SetSelectedPosition(p3New); this->GetRenderWindowPart()->RequestUpdate(); } } } void QmitkTensorReconstructionView::Advanced1CheckboxClicked() { bool check = m_Controls-> m_Advanced1->isChecked(); m_Controls->frame->setVisible(check); } void QmitkTensorReconstructionView::Activated() { } void QmitkTensorReconstructionView::Deactivated() { // Get all current nodes mitk::DataStorage::SetOfObjects::ConstPointer objects = this->GetDataStorage()->GetAll(); mitk::DataStorage::SetOfObjects::const_iterator itemiter( objects->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( objects->end() ); while ( itemiter != itemiterend ) // for all items { mitk::DataNode::Pointer node = *itemiter; if (node.IsNull()) continue; // only look at interesting types if( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(dynamic_cast(node->GetData()))) { if (this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", *itemiter)) { node = this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", *itemiter); this->GetDataStorage()->Remove(node); } } itemiter++; } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkTensorReconstructionView::Visible() { } void QmitkTensorReconstructionView::Hidden() { } void QmitkTensorReconstructionView::ResidualCalculation() { // Extract dwi and dti from current selection // In case of multiple selections, take the first one, since taking all combinations is not meaningful mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); mitk::Image::Pointer diffImage = mitk::Image::New(); TensorImageType::Pointer tensorImage; std::string nodename; if(m_Controls->m_ResBox1->GetSelectedNode()) { diffImage = static_cast(m_Controls->m_ResBox1->GetSelectedNode()->GetData()); } else return; if(m_Controls->m_ResBox2->GetSelectedNode().IsNotNull()) { mitk::TensorImage* mitkVol; mitkVol = static_cast(m_Controls->m_ResBox2->GetSelectedNode()->GetData()); mitk::CastToItkImage(mitkVol, tensorImage); m_Controls->m_ResBox2->GetSelectedNode()->GetStringProperty("name", nodename); } else return; typedef itk::TensorImageToDiffusionImageFilter FilterType; GradientDirectionContainerType::ConstPointer gradients = mitk::DiffusionPropertyHelper::GetGradientContainer(diffImage); // Find the min and the max values from a baseline image mitk::ImageStatisticsHolder *stats = diffImage->GetStatistics(); //Initialize filter that calculates the modeled diffusion weighted signals FilterType::Pointer filter = FilterType::New(); filter->SetInput( tensorImage ); filter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(diffImage)); filter->SetGradientList(gradients); filter->SetMin(stats->GetScalarValueMin()); filter->SetMax(stats->GetScalarValueMax()); filter->Update(); mitk::Image::Pointer image = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::CopyProperties(diffImage, image, true); mitk::DiffusionPropertyHelper::SetGradientContainer(image, gradients); mitk::DiffusionPropertyHelper::InitializeImage( image ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); mitk::ImageVtkMapper2D::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_Estimated DWI"); node->SetName(newname.toLatin1()); GetDataStorage()->Add(node, m_Controls->m_ResBox2->GetSelectedNode()); BValueMapType map = mitk::DiffusionPropertyHelper::GetBValueMap(image); std::vector< unsigned int > b0Indices = map[0]; typedef itk::ResidualImageFilter ResidualImageFilterType; ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(diffImage, itkVectorImagePointer); ITKDiffusionImageType::Pointer itkSecondVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(image, itkSecondVectorImagePointer); ResidualImageFilterType::Pointer residualFilter = ResidualImageFilterType::New(); residualFilter->SetInput( itkVectorImagePointer ); residualFilter->SetSecondDiffusionImage( itkSecondVectorImagePointer ); residualFilter->SetGradients(gradients); residualFilter->SetB0Index(b0Indices[0]); residualFilter->SetB0Threshold(30); residualFilter->Update(); itk::Image::Pointer residualImage = itk::Image::New(); residualImage = residualFilter->GetOutput(); mitk::Image::Pointer mitkResImg = mitk::Image::New(); mitk::CastToMitkImage(residualImage, mitkResImg); stats = mitkResImg->GetStatistics(); float min = stats->GetScalarValueMin(); float max = stats->GetScalarValueMax(); mitk::LookupTableProperty::Pointer lutProp = mitk::LookupTableProperty::New(); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(min, max); // If you don't want to use the whole color range, you can use // SetValueRange, SetHueRange, and SetSaturationRange lookupTable->Build(); vtkSmartPointer reversedlookupTable = vtkSmartPointer::New(); reversedlookupTable->SetTableRange(min+1, max); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookupTable->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } lut->SetVtkLookupTable(lookupTable); lutProp->SetLookupTable(lut); // Create lookuptable mitk::DataNode::Pointer resNode=mitk::DataNode::New(); resNode->SetData( mitkResImg ); resNode->SetName("Residuals"); resNode->SetProperty("LookupTable", lutProp); bool b; resNode->GetBoolProperty("use color", b); resNode->SetBoolProperty("use color", false); GetDataStorage()->Add(resNode, m_Controls->m_ResBox2->GetSelectedNode()); this->GetRenderWindowPart()->RequestUpdate(); // Draw Graph std::vector means = residualFilter->GetMeans(); std::vector q1s = residualFilter->GetQ1(); std::vector q3s = residualFilter->GetQ3(); std::vector percentagesOfOUtliers = residualFilter->GetPercentagesOfOutliers(); m_Controls->m_ResidualAnalysis->SetMeans(means); m_Controls->m_ResidualAnalysis->SetQ1(q1s); m_Controls->m_ResidualAnalysis->SetQ3(q3s); m_Controls->m_ResidualAnalysis->SetPercentagesOfOutliers(percentagesOfOUtliers); if(m_Controls->m_PercentagesOfOutliers->isChecked()) { m_Controls->m_ResidualAnalysis->DrawPercentagesOfOutliers(); } else { m_Controls->m_ResidualAnalysis->DrawMeans(); } // Draw Graph for volumes per slice in the QGraphicsView std::vector< std::vector > outliersPerSlice = residualFilter->GetOutliersPerSlice(); int xSize = outliersPerSlice.size(); if(xSize == 0) { return; } int ySize = outliersPerSlice[0].size(); // Find maximum in outliersPerSlice double maxOutlier= 0.0; for(int i=0; imaxOutlier) { maxOutlier = outliersPerSlice[i][j]; } } } // Create some QImage QImage qImage(xSize, ySize, QImage::Format_RGB32); QImage legend(1, 256, QImage::Format_RGB32); QRgb value; vtkSmartPointer lookup = vtkSmartPointer::New(); lookup->SetTableRange(0.0, maxOutlier); lookup->Build(); reversedlookupTable->SetTableRange(0, maxOutlier); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookup->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } // Fill qImage for(int i=0; iMapValue(out); + const unsigned char *_rgba = lookup->MapValue(out); int r, g, b; r = _rgba[0]; g = _rgba[1]; b = _rgba[2]; value = qRgb(r, g, b); qImage.setPixel(i,j,value); } } for(int i=0; i<256; i++) { double* rgba = lookup->GetTableValue(i); int r, g, b; r = rgba[0]*255; g = rgba[1]*255; b = rgba[2]*255; value = qRgb(r, g, b); legend.setPixel(0,255-i,value); } QString upper = QString::number(maxOutlier, 'g', 3); upper.append(" %"); QString lower = QString::number(0.0); lower.append(" %"); m_Controls->m_UpperLabel->setText(upper); m_Controls->m_LowerLabel->setText(lower); QPixmap pixmap(QPixmap::fromImage(qImage)); QGraphicsPixmapItem *item = new QGraphicsPixmapItem(pixmap); item->setTransform(QTransform::fromScale(10.0, 3.0), true); QPixmap pixmap2(QPixmap::fromImage(legend)); QGraphicsPixmapItem *item2 = new QGraphicsPixmapItem(pixmap2); item2->setTransform(QTransform::fromScale(20.0, 1.0), true); m_Controls->m_PerSliceView->SetResidualPixmapItem(item); QGraphicsScene* scene = new QGraphicsScene; QGraphicsScene* scene2 = new QGraphicsScene; scene->addItem(item); scene2->addItem(item2); m_Controls->m_PerSliceView->setScene(scene); m_Controls->m_LegendView->setScene(scene2); m_Controls->m_PerSliceView->show(); m_Controls->m_PerSliceView->repaint(); m_Controls->m_LegendView->setHorizontalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->setVerticalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->show(); m_Controls->m_LegendView->repaint(); } void QmitkTensorReconstructionView::Reconstruct() { int method = m_Controls->m_ReconctructionMethodBox->currentIndex(); switch (method) { case 0: ItkTensorReconstruction(); break; case 1: TensorReconstructionWithCorr(); break; default: ItkTensorReconstruction(); } } void QmitkTensorReconstructionView::TensorReconstructionWithCorr() { try { if ( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ) // for all items { mitk::Image* vols = static_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); std::string nodename = m_Controls->m_DwiBox->GetSelectedNode()->GetName(); typedef itk::TensorReconstructionWithEigenvalueCorrectionFilter< DiffusionPixelType, TTensorPixelType > ReconstructionFilter; float b0Threshold = m_Controls->m_TensorReconstructionThreshold->value(); mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientContainerCopy = mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New(); for(auto it = mitk::DiffusionPropertyHelper::GetGradientContainer(vols)->Begin(); it != mitk::DiffusionPropertyHelper::GetGradientContainer(vols)->End(); it++) gradientContainerCopy->push_back(it.Value()); ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); ReconstructionFilter::Pointer reconFilter = ReconstructionFilter::New(); reconFilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(vols)); reconFilter->SetGradientImage( mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::ConstPointer(gradientContainerCopy), itkVectorImagePointer); reconFilter->SetB0Threshold(b0Threshold); reconFilter->Update(); typedef mitk::TensorImage::ItkTensorImageType TensorImageType; TensorImageType::Pointer outputTensorImg = reconFilter->GetOutput(); typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(outputTensorImg, outputTensorImg->GetRequestedRegion()); tensorIt.GoToBegin(); int negatives = 0; while(!tensorIt.IsAtEnd()) { typedef mitk::TensorImage::PixelType TensorType; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iInitializeByItk( outputTensorImg.GetPointer() ); image->SetVolume( outputTensorImg->GetBufferPointer() ); mitk::DiffusionPropertyHelper::CopyDICOMProperties(vols, image); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetName(nodename+"_EigenvalueCorrected_DT"); GetDataStorage()->Add(node, m_Controls->m_DwiBox->GetSelectedNode()); } this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); } } void QmitkTensorReconstructionView::ItkTensorReconstruction() { try { if ( m_Controls->m_DwiBox->GetSelectedNode().IsNotNull() ) // for all items { mitk::Image* vols = static_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); std::string nodename = m_Controls->m_DwiBox->GetSelectedNode()->GetName(); typedef itk::DiffusionTensor3DReconstructionImageFilter TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = mitk::DiffusionPropertyHelper::GetGradientContainer(vols)->Begin(); it != mitk::DiffusionPropertyHelper::GetGradientContainer(vols)->End(); it++) { gradientContainerCopy->push_back(it.Value()); } ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); tensorReconstructionFilter->SetBValue(mitk::DiffusionPropertyHelper::GetReferenceBValue(vols)); tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, itkVectorImagePointer ); tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreshold->value() ); tensorReconstructionFilter->Update(); // TENSORS TO DATATREE mitk::TensorImage::Pointer image = mitk::TensorImage::New(); typedef mitk::TensorImage::ItkTensorImageType TensorImageType; TensorImageType::Pointer tensorImage; tensorImage = tensorReconstructionFilter->GetOutput(); // Check the tensor for negative eigenvalues if(m_Controls->m_CheckNegativeEigenvalues->isChecked()) { typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetRequestedRegion()); tensorIt.GoToBegin(); while(!tensorIt.IsAtEnd()) { typedef mitk::TensorImage::PixelType TensorType; //typedef itk::Tensor TensorType2; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iSetDirection( itkVectorImagePointer->GetDirection() ); image->InitializeByItk( tensorImage.GetPointer() ); image->SetVolume( tensorReconstructionFilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); mitk::DiffusionPropertyHelper::CopyDICOMProperties(vols, image); node->SetData( image ); node->SetName(nodename+"_LinearLeastSquares_DT"); GetDataStorage()->Add(node, m_Controls->m_DwiBox->GetSelectedNode()); } this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return; } } void QmitkTensorReconstructionView::TensorsToDWI() { DoTensorsToDWI(); } void QmitkTensorReconstructionView::TensorsToOdf() { if (m_Controls->m_OdfBox->GetSelectedNode().IsNotNull()) { mitk::DataNode::Pointer tensorImageNode = m_Controls->m_OdfBox->GetSelectedNode(); typedef mitk::TensorImage::ScalarPixelType TTensorPixelType; typedef mitk::TensorImage::ItkTensorImageType TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(dynamic_cast(tensorImageNode->GetData()), itkvol); typedef itk::TensorImageToOdfImageFilter< TTensorPixelType, TTensorPixelType > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->Update(); typedef itk::Vector OutputPixelType; typedef itk::Image OutputImageType; mitk::OdfImage::Pointer image = mitk::OdfImage::New(); OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); mitk::DiffusionPropertyHelper::CopyDICOMProperties(tensorImageNode->GetData(), image); node->SetData( image ); node->SetName(tensorImageNode->GetName()+"_Odf"); GetDataStorage()->Add(node, tensorImageNode); } } void QmitkTensorReconstructionView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { UpdateGui(); } void QmitkTensorReconstructionView::UpdateGui() { m_Controls->m_StartReconstruction->setEnabled(m_Controls->m_DwiBox->GetSelectedNode().IsNotNull()); m_Controls->m_TensorsToDWIButton->setEnabled(m_Controls->m_DtiBox->GetSelectedNode().IsNotNull()); m_Controls->m_TensorsToOdfButton->setEnabled(m_Controls->m_OdfBox->GetSelectedNode().IsNotNull()); m_Controls->m_ResidualButton->setEnabled(m_Controls->m_ResBox1->GetSelectedNode().IsNotNull() && m_Controls->m_ResBox2->GetSelectedNode().IsNotNull()); m_Controls->m_PercentagesOfOutliers->setEnabled(m_Controls->m_ResBox1->GetSelectedNode().IsNotNull() && m_Controls->m_ResBox2->GetSelectedNode().IsNotNull()); m_Controls->m_PerSliceView->setEnabled(m_Controls->m_ResBox1->GetSelectedNode().IsNotNull() && m_Controls->m_ResBox2->GetSelectedNode().IsNotNull()); } template itk::VectorContainer >::Pointer QmitkTensorReconstructionView::MakeGradientList() { itk::VectorContainer >::Pointer retval = itk::VectorContainer >::New(); vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); for(int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval->push_back(v); } // Add 0 vector for B0 vnl_vector_fixed v; v.fill(0.0); retval->push_back(v); return retval; } void QmitkTensorReconstructionView::DoTensorsToDWI() { try { if (m_Controls->m_DtiBox->GetSelectedNode().IsNotNull()) { std::string nodename = m_Controls->m_DtiBox->GetSelectedNode()->GetName(); mitk::TensorImage* vol = static_cast(m_Controls->m_DtiBox->GetSelectedNode()->GetData()); typedef mitk::TensorImage::ScalarPixelType TTensorPixelType; typedef mitk::TensorImage::ItkTensorImageType TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(vol, itkvol); typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; FilterType::GradientListPointerType gradientList; switch(m_Controls->m_TensorsToDWINumDirsSelect->currentIndex()) { case 0: gradientList = MakeGradientList<12>(); break; case 1: gradientList = MakeGradientList<42>(); break; case 2: gradientList = MakeGradientList<92>(); break; case 3: gradientList = MakeGradientList<162>(); break; case 4: gradientList = MakeGradientList<252>(); break; case 5: gradientList = MakeGradientList<362>(); break; case 6: gradientList = MakeGradientList<492>(); break; case 7: gradientList = MakeGradientList<642>(); break; case 8: gradientList = MakeGradientList<812>(); break; case 9: gradientList = MakeGradientList<1002>(); break; default: gradientList = MakeGradientList<92>(); } double bVal = m_Controls->m_TensorsToDWIBValueEdit->text().toDouble(); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->SetBValue(bVal); filter->SetGradientList(gradientList); filter->Update(); mitk::Image::Pointer image = mitk::GrabItkImageMemory( filter->GetOutput() ); mitk::DiffusionPropertyHelper::SetGradientContainer(image, gradientList); mitk::DiffusionPropertyHelper::SetReferenceBValue(image, bVal); mitk::DiffusionPropertyHelper::InitializeImage( image ); mitk::DataNode::Pointer node=mitk::DataNode::New(); mitk::DiffusionPropertyHelper::CopyDICOMProperties(vol, image); node->SetData( image ); mitk::ImageVtkMapper2D::SetDefaultProperties(node); node->SetName(nodename+"_DWI"); GetDataStorage()->Add(node, m_Controls->m_DtiBox->GetSelectedNode()); } this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "DWI estimation failed:", ex.GetDescription()); return ; } } void QmitkTensorReconstructionView::PreviewThreshold(int threshold) { if (m_Controls->m_DwiBox->GetSelectedNode().IsNotNull()) { mitk::Image* vols = static_cast(m_Controls->m_DwiBox->GetSelectedNode()->GetData()); ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); // Extract b0 image typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filterB0 = FilterType::New(); filterB0->SetInput(itkVectorImagePointer); filterB0->SetDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(vols)); filterB0->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); typedef itk::Image ImageType; typedef itk::Image SegmentationType; typedef itk::BinaryThresholdImageFilter ThresholdFilterType; // apply threshold ThresholdFilterType::Pointer filterThreshold = ThresholdFilterType::New(); filterThreshold->SetInput(filterB0->GetOutput()); filterThreshold->SetLowerThreshold(threshold); filterThreshold->SetInsideValue(0); filterThreshold->SetOutsideValue(1); // mark cut off values red filterThreshold->Update(); mitkImage->InitializeByItk( filterThreshold->GetOutput() ); mitkImage->SetVolume( filterThreshold->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node; if (this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", m_Controls->m_DwiBox->GetSelectedNode())) { node = this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", m_Controls->m_DwiBox->GetSelectedNode()); } else { // create a new node, to show thresholded values node = mitk::DataNode::New(); GetDataStorage()->Add( node, m_Controls->m_DwiBox->GetSelectedNode() ); node->SetProperty( "name", mitk::StringProperty::New("ThresholdOverlay")); node->SetBoolProperty("helper object", true); } node->SetData( mitkImage ); mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp index 227276e..c94cfdb 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp @@ -1,606 +1,601 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkGibbsTrackingView.h" // Qt #include #include #include #include // MITK #include #include #include #include #include #include #include #include #include #include // ITK #include #include #include // MISC -#include +#include QmitkTrackingWorker::QmitkTrackingWorker(QmitkGibbsTrackingView* view) : m_View(view) { } void QmitkTrackingWorker::run() { m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); m_View->m_GlobalTracker->SetOdfImage(m_View->m_ItkOdfImage); m_View->m_GlobalTracker->SetTensorImage(m_View->m_ItkTensorImage); m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); m_View->m_GlobalTracker->SetStartTemperature((float)m_View->m_Controls->m_StartTempSlider->value()/100); m_View->m_GlobalTracker->SetEndTemperature((float)m_View->m_Controls->m_EndTempSlider->value()/10000); m_View->m_GlobalTracker->SetIterations(m_View->m_Controls->m_IterationsBox->text().toDouble()); m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); m_View->m_GlobalTracker->SetMinFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); m_View->m_GlobalTracker->SetCurvatureThreshold(cos((float)m_View->m_Controls->m_CurvatureThresholdSlider->value()*itk::Math::pi/180)); m_View->m_GlobalTracker->SetRandomSeed(m_View->m_Controls->m_RandomSeedSlider->value()); try{ m_View->m_GlobalTracker->Update(); } catch( mitk::Exception& e ) { MITK_ERROR << "Internal error occured: " << e.what() << "\nAborting"; } m_View->m_TrackingThread.quit(); } const std::string QmitkGibbsTrackingView::VIEW_ID = "org.mitk.views.gibbstracking"; QmitkGibbsTrackingView::QmitkGibbsTrackingView() : QmitkAbstractView() , m_Controls( 0 ) , m_TrackingNode(nullptr) , m_FiberBundle(nullptr) , m_MaskImage(nullptr) , m_TensorImage(nullptr) , m_OdfImage(nullptr) , m_ItkOdfImage(nullptr) , m_ItkTensorImage(nullptr) , m_ImageNode(nullptr) , m_MaskImageNode(nullptr) , m_FiberBundleNode(nullptr) , m_ThreadIsRunning(false) , m_ElapsedTime(0) , m_GlobalTracker(nullptr) , m_TrackingWorker(this) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } QmitkGibbsTrackingView::~QmitkGibbsTrackingView() { if (m_GlobalTracker.IsNull()) return; m_GlobalTracker->SetAbortTracking(true); m_TrackingThread.wait(); } // update tracking status and generate fiber bundle void QmitkGibbsTrackingView::TimerUpdate() { UpdateTrackingStatus(); GenerateFiberBundle(); } // tell global tractography filter to stop after current step void QmitkGibbsTrackingView::StopGibbsTracking() { if (m_GlobalTracker.IsNull()) return; m_GlobalTracker->SetAbortTracking(true); m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); m_TrackingNode = nullptr; } // update gui elements and generate fiber bundle after tracking is finished void QmitkGibbsTrackingView::AfterThread() { m_ThreadIsRunning = false; m_TrackingTimer->stop(); UpdateGUI(); if( !m_GlobalTracker->GetIsInValidState() ) { QMessageBox::critical( nullptr, "Gibbs Tracking", "An internal error occured. Tracking aborted.\n Please check the log for details." ); m_FiberBundleNode = nullptr; return; } UpdateTrackingStatus(); if(m_Controls->m_ParticleWeightSlider->value()==0) { m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); } if(m_Controls->m_ParticleWidthSlider->value()==0) { m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); } if(m_Controls->m_ParticleLengthSlider->value()==0) { m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); } GenerateFiberBundle(); m_FiberBundleNode = 0; m_GlobalTracker = 0; // images not needed anymore ( relevant only for computation ) // we need to release them to remove the memory access block created through CastToItk<> calls this->m_ItkOdfImage = 0; this->m_ItkTensorImage = 0; } // start tracking timer and update gui elements before tracking is started void QmitkGibbsTrackingView::BeforeThread() { m_ThreadIsRunning = true; m_TrackingTime = QTime::currentTime(); m_ElapsedTime = 0; m_TrackingTimer->start(1000); UpdateGUI(); } // setup gui elements and signal/slot connections void QmitkGibbsTrackingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkGibbsTrackingViewControls; m_Controls->setupUi( parent ); AdvancedSettings(); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); connect( m_Controls->m_CurvatureThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(SetCurvatureThreshold(int)) ); connect( m_Controls->m_RandomSeedSlider, SIGNAL(valueChanged(int)), this, SLOT(SetRandomSeed(int)) ); connect( m_Controls->m_OutputFileButton, SIGNAL(clicked()), this, SLOT(SetOutputFile()) ); m_Controls->m_InputImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isOdfImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTensorImagePredicate = mitk::TNodePredicateDataType::New(); m_Controls->m_InputImageBox->SetPredicate( mitk::NodePredicateOr::New(isOdfImagePredicate, isTensorImagePredicate) ); m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateDimension::Pointer is3D = mitk::NodePredicateDimension::New(3); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, mitk::NodePredicateAnd::New(isImagePredicate, is3D)) ); connect( (QObject*)(m_Controls->m_MaskImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGUI())); connect( (QObject*)(m_Controls->m_InputImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGUI())); } } void QmitkGibbsTrackingView::SetFocus() { m_Controls->m_TrackingStart->setFocus(); } void QmitkGibbsTrackingView::SetInExBalance(int value) { m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); } void QmitkGibbsTrackingView::SetFiberLength(int value) { m_Controls->m_FiberLengthLabel->setText(QString::number(value)+"mm"); } void QmitkGibbsTrackingView::SetRandomSeed(int value) { if (value>=0) m_Controls->m_RandomSeedLabel->setText(QString::number(value)); else m_Controls->m_RandomSeedLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleWeight(int value) { if (value>0) m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); else m_Controls->m_ParticleWeightLabel->setText("auto"); } void QmitkGibbsTrackingView::SetStartTemp(int value) { m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); } void QmitkGibbsTrackingView::SetEndTemp(int value) { m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); } void QmitkGibbsTrackingView::SetParticleWidth(int value) { if (value>0) m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleWidthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleLength(int value) { if (value>0) m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleLengthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetCurvatureThreshold(int value) { m_Controls->m_CurvatureThresholdLabel->setText(QString::number(value)+"°"); } // called if datamanager selection changes void QmitkGibbsTrackingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { UpdateGUI(); } void QmitkGibbsTrackingView::NodeRemoved(const mitk::DataNode * node) { if (m_ThreadIsRunning) { if (node==m_TrackingNode.GetPointer()) { StopGibbsTracking(); } } } // update gui elements displaying trackings status void QmitkGibbsTrackingView::UpdateTrackingStatus() { if (m_GlobalTracker.IsNull()) return; m_ElapsedTime += m_TrackingTime.elapsed()/1000; m_TrackingTime.restart(); unsigned long hours = m_ElapsedTime/3600; unsigned long minutes = (m_ElapsedTime%3600)/60; unsigned long seconds = m_ElapsedTime%60; m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); m_Controls->m_CurrentStepLabel->setText( QString::number(100*m_GlobalTracker->GetCurrentIteration()/m_GlobalTracker->GetIterations())+"%" ); m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); } // update gui elements (enable/disable elements and set tooltips) void QmitkGibbsTrackingView::UpdateGUI() { if (!m_ThreadIsRunning && m_Controls->m_InputImageBox->GetSelectedNode().IsNotNull()) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(true); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_IterationsBox->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); m_Controls->m_TrackingStop->setToolTip(""); m_Controls->m_MaskImageBox->setEnabled(true); m_Controls->m_InputImageBox->setEnabled(true); } else if (!m_ThreadIsRunning) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_IterationsBox->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("No ODF image selected."); m_Controls->m_TrackingStop->setToolTip(""); m_Controls->m_MaskImageBox->setEnabled(true); m_Controls->m_InputImageBox->setEnabled(true); } else { m_Controls->m_TrackingStop->setEnabled(true); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(false); m_Controls->m_IterationsBox->setEnabled(false); m_Controls->m_AdvancedFrame->setEnabled(false); m_Controls->m_AdvancedFrame->setVisible(false); m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); m_Controls->m_MaskImageBox->setEnabled(false); m_Controls->m_InputImageBox->setEnabled(false); } } // show/hide advanced settings frame void QmitkGibbsTrackingView::AdvancedSettings() { m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); } // check for mask and odf and start tracking thread void QmitkGibbsTrackingView::StartGibbsTracking() { if(m_ThreadIsRunning) { MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; return; } m_GlobalTracker = nullptr; if (m_Controls->m_InputImageBox->GetSelectedNode().IsNull()) { QMessageBox::information( nullptr, "Warning", "Please load and select a Odf image before starting image processing."); return; } m_ImageNode = m_Controls->m_InputImageBox->GetSelectedNode(); if (dynamic_cast(m_ImageNode->GetData())) m_OdfImage = dynamic_cast(m_ImageNode->GetData()); else if (dynamic_cast(m_ImageNode->GetData())) m_TensorImage = dynamic_cast(m_ImageNode->GetData()); if (m_OdfImage.IsNull() && m_TensorImage.IsNull()) return; // cast odf to itk m_TrackingNode = m_ImageNode; m_ItkTensorImage = nullptr; m_ItkOdfImage = nullptr; m_MaskImage = nullptr; if (m_OdfImage.IsNotNull()) { m_ItkOdfImage = ItkOdfImgType::New(); mitk::CastToItkImage(m_OdfImage, m_ItkOdfImage); } else { m_ItkTensorImage = ItkTensorImage::New(); mitk::CastToItkImage(m_TensorImage, m_ItkTensorImage); } // mask image found? // catch exceptions thrown by the itkAccess macros try{ if(m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) { m_MaskImageNode = m_Controls->m_MaskImageBox->GetSelectedNode(); if (dynamic_cast(m_MaskImageNode->GetData())) mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); } } catch(...){}; // start worker thread m_TrackingThread.start(QThread::LowestPriority); } // generate mitkFiberBundle from tracking filter output void QmitkGibbsTrackingView::GenerateFiberBundle() { if (m_GlobalTracker.IsNull() || (!(m_Controls->m_VisualizationCheckbox->isChecked() || m_Controls->m_VisualizeOnceButton->isChecked()) && m_ThreadIsRunning)) return; if (m_Controls->m_VisualizeOnceButton->isChecked()) m_Controls->m_VisualizeOnceButton->setChecked(false); vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); if ( m_GlobalTracker->GetNumAcceptedFibers()==0 ) return; m_FiberBundle = mitk::FiberBundle::New(fiberBundle); m_FiberBundle->SetTrackVisHeader(dynamic_cast(m_ImageNode->GetData())->GetGeometry()); mitk::DiffusionPropertyHelper::CopyDICOMProperties(m_ImageNode->GetData(), m_FiberBundle); if (m_FiberBundleNode.IsNotNull()){ GetDataStorage()->Remove(m_FiberBundleNode); m_FiberBundleNode = 0; } m_GlobalTracker->SetDicomProperties(m_FiberBundle); m_FiberBundleNode = mitk::DataNode::New(); m_FiberBundleNode->SetData(m_FiberBundle); QString name("FiberBundle_"); name += m_ImageNode->GetName().c_str(); name += "_Gibbs"; m_FiberBundleNode->SetName(name.toStdString()); m_FiberBundleNode->SetVisibility(true); if (!m_OutputFileName.isEmpty() && !m_ThreadIsRunning) { try { mitk::IOUtil::Save(m_FiberBundle.GetPointer(),m_OutputFileName.toStdString()); QMessageBox::information(nullptr, "Fiber bundle saved to", m_OutputFileName); } catch (itk::ExceptionObject &ex) { QMessageBox::information(nullptr, "Fiber bundle could not be saved", QString("%1\n%2\n%3\n%4\n%5\n%6").arg(ex.GetNameOfClass()).arg(ex.GetFile()).arg(ex.GetLine()).arg(ex.GetLocation()).arg(ex.what()).arg(ex.GetDescription())); } } GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); } void QmitkGibbsTrackingView::SetOutputFile() { // SELECT FOLDER DIALOG m_OutputFileName = QFileDialog::getSaveFileName(0, tr("Set file name"), QDir::currentPath()+"/FiberBundle.fib", tr("Fiber Bundle (*.fib)") ); if (m_OutputFileName.isEmpty()) m_Controls->m_OutputFileLabel->setText("N/A"); else m_Controls->m_OutputFileLabel->setText(m_OutputFileName); } // save current tracking paramters as xml file (.gtp) void QmitkGibbsTrackingView::SaveTrackingParameters() { - TiXmlDocument documentXML; - TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); - documentXML.LinkEndChild( declXML ); - - TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); - mainXML->SetAttribute("file_version", "0.1"); - documentXML.LinkEndChild(mainXML); - - TiXmlElement* paramXML = new TiXmlElement("parameter_set"); - paramXML->SetAttribute("iterations", m_Controls->m_IterationsBox->text().toStdString()); - paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); - paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); - paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); - paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); - paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); - paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); - paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); - paramXML->SetAttribute("curvature_threshold", QString::number(m_Controls->m_CurvatureThresholdSlider->value()).toStdString()); - mainXML->LinkEndChild(paramXML); + tinyxml2::XMLDocument documentXML; + //TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); + //documentXML.LinkEndChild( declXML ); + + auto mainXML = documentXML.NewElement("global_tracking_parameter_file"); + mainXML->SetAttribute("file_version", 0.1); + documentXML.InsertFirstChild(mainXML); + + auto paramXML = documentXML.NewElement("parameter_set"); + paramXML->SetAttribute("iterations", m_Controls->m_IterationsBox->text().toStdString().c_str()); + paramXML->SetAttribute("particle_length", (float)m_Controls->m_ParticleLengthSlider->value()/10); + paramXML->SetAttribute("particle_width", (float)m_Controls->m_ParticleWidthSlider->value()/10); + paramXML->SetAttribute("particle_weight", (float)m_Controls->m_ParticleWeightSlider->value()/10000); + paramXML->SetAttribute("temp_start", (float)m_Controls->m_StartTempSlider->value()/100); + paramXML->SetAttribute("temp_end", (float)m_Controls->m_EndTempSlider->value()/10000); + paramXML->SetAttribute("inexbalance", (float)m_Controls->m_InExBalanceSlider->value()/10); + paramXML->SetAttribute("fiber_length", m_Controls->m_FiberLengthSlider->value()); + paramXML->SetAttribute("curvature_threshold", m_Controls->m_CurvatureThresholdSlider->value()); + mainXML->InsertEndChild(paramXML); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), QDir::currentPath()+"/param.gtp", tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; if(!filename.endsWith(".gtp")) filename += ".gtp"; - documentXML.SaveFile( filename.toStdString() ); + documentXML.SaveFile( filename.toStdString().c_str() ); } // load current tracking paramters from xml file (.gtp) void QmitkGibbsTrackingView::LoadTrackingParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; - TiXmlDocument doc( filename.toStdString() ); - doc.LoadFile(); + tinyxml2::XMLDocument doc; + doc.LoadFile(filename.toStdString().c_str()); - TiXmlHandle hDoc(&doc); - TiXmlElement* pElem; - TiXmlHandle hRoot(0); - - pElem = hDoc.FirstChildElement().Element(); - hRoot = TiXmlHandle(pElem); - pElem = hRoot.FirstChildElement("parameter_set").Element(); + tinyxml2::XMLNode * hRoot = doc.FirstChild(); + auto pElem = hRoot->FirstChildElement("parameter_set"); QString iterations(pElem->Attribute("iterations")); m_Controls->m_IterationsBox->setText(iterations); QString particleLength(pElem->Attribute("particle_length")); float pLength = particleLength.toFloat(); QString particleWidth(pElem->Attribute("particle_width")); float pWidth = particleWidth.toFloat(); if (pLength==0) m_Controls->m_ParticleLengthLabel->setText("auto"); else m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); if (pWidth==0) m_Controls->m_ParticleWidthLabel->setText("auto"); else m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); m_Controls->m_ParticleLengthSlider->setValue(pLength*10); QString partWeight(pElem->Attribute("particle_weight")); m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); m_Controls->m_ParticleWeightLabel->setText(partWeight); QString startTemp(pElem->Attribute("temp_start")); m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); m_Controls->m_StartTempLabel->setText(startTemp); QString endTemp(pElem->Attribute("temp_end")); m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); m_Controls->m_EndTempLabel->setText(endTemp); QString inExBalance(pElem->Attribute("inexbalance")); m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); m_Controls->m_InExBalanceLabel->setText(inExBalance); QString fiberLength(pElem->Attribute("fiber_length")); m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); m_Controls->m_FiberLengthLabel->setText(fiberLength+"mm"); QString curvThres(pElem->Attribute("curvature_threshold")); m_Controls->m_CurvatureThresholdSlider->setValue(curvThres.toInt()); m_Controls->m_CurvatureThresholdLabel->setText(curvThres+"°"); }