diff --git a/Applications/Diffusion/CMakeLists.txt b/Applications/Diffusion/CMakeLists.txt index 903b7d2193..fc47284c90 100644 --- a/Applications/Diffusion/CMakeLists.txt +++ b/Applications/Diffusion/CMakeLists.txt @@ -1,95 +1,96 @@ project(MitkDiffusion) set(DIFFUSIONAPP_NAME MitkDiffusion) set(_app_options) if(MITK_SHOW_CONSOLE_WINDOW) list(APPEND _app_options SHOW_CONSOLE) endif() # Create a cache entry for the provisioning file which is used to export # the file name in the MITKConfig.cmake file. This will keep external projects # which rely on this file happy. set(DIFFUSIONIMAGINGAPP_PROVISIONING_FILE "${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${DIFFUSIONAPP_NAME}.provisioning" CACHE INTERNAL "${DIFFUSIONAPP_NAME} provisioning file" FORCE) # should be identical to the list in /CMake/mitkBuildConfigurationmitkDiffusion.cmake # remember to set plugins which should be automatically toggled in target_libraries.cmake set(_plugins org.commontk.configadmin org.commontk.eventadmin org.blueberry.core.runtime org.blueberry.core.expressions org.blueberry.core.commands org.blueberry.ui.qt org.blueberry.ui.qt.log org.blueberry.ui.qt.help org.mitk.core.services org.mitk.gui.common org.mitk.planarfigure org.mitk.core.ext org.mitk.gui.qt.application org.mitk.gui.qt.ext org.mitk.gui.qt.diffusionimagingapp org.mitk.gui.qt.common org.mitk.gui.qt.stdmultiwidgeteditor org.mitk.gui.qt.common.legacy org.mitk.gui.qt.datamanager org.mitk.gui.qt.measurementtoolbox org.mitk.gui.qt.segmentation org.mitk.gui.qt.volumevisualization org.mitk.gui.qt.diffusionimaging org.mitk.gui.qt.diffusionimaging.connectomics org.mitk.gui.qt.diffusionimaging.fiberfox org.mitk.gui.qt.diffusionimaging.fiberprocessing org.mitk.gui.qt.diffusionimaging.ivim org.mitk.gui.qt.diffusionimaging.odfpeaks org.mitk.gui.qt.diffusionimaging.preprocessing org.mitk.gui.qt.diffusionimaging.reconstruction org.mitk.gui.qt.diffusionimaging.tractography org.mitk.gui.qt.matchpoint.algorithm.browser org.mitk.gui.qt.matchpoint.algorithm.control + org.mitk.gui.qt.matchpoint.mapper org.mitk.gui.qt.imagenavigator org.mitk.gui.qt.moviemaker org.mitk.gui.qt.basicimageprocessing org.mitk.gui.qt.properties org.mitk.gui.qt.viewnavigator ) # Plug-ins listed below will not be # - added as a build-time dependency to the executable # - listed in the provisioning file for the executable # - installed if they are external plug-ins set(_exclude_plugins org.blueberry.test org.blueberry.uitest org.mitk.gui.qt.coreapplication org.mitk.gui.qt.extapplication ) set(_src_files MitkDiffusion.cpp ) qt5_add_resources(_src_files splashscreen.qrc) mitkFunctionCreateBlueBerryApplication( NAME ${DIFFUSIONAPP_NAME} DESCRIPTION "MITK Diffusion" PLUGINS ${_plugins} EXCLUDE_PLUGINS ${_exclude_plugins} SOURCES ${_src_files} ${_app_options} ) mitk_use_modules(TARGET ${DIFFUSIONAPP_NAME} MODULES MitkAppUtil) # Add meta dependencies (e.g. on auto-load modules from depending modules) if(TARGET ${CMAKE_PROJECT_NAME}-autoload) add_dependencies(${DIFFUSIONAPP_NAME} ${CMAKE_PROJECT_NAME}-autoload) endif() # Add a build time dependency to legacy BlueBerry bundles. if(MITK_MODULES_ENABLED_PLUGINS) add_dependencies(${DIFFUSIONAPP_NAME} ${MITK_MODULES_ENABLED_PLUGINS}) endif() diff --git a/Applications/Diffusion/target_libraries.cmake b/Applications/Diffusion/target_libraries.cmake index a6496917be..347d5fbc04 100644 --- a/Applications/Diffusion/target_libraries.cmake +++ b/Applications/Diffusion/target_libraries.cmake @@ -1,29 +1,32 @@ # A list of plug-in targets which should be automatically enabled -# (or be available in external projects) for this application. +# (or be available in external projects) for this application set(target_libraries org_blueberry_ui_qt org_blueberry_ui_qt_help org_mitk_planarfigure org_mitk_gui_qt_diffusionimagingapp org_mitk_gui_qt_common_legacy org_mitk_gui_qt_ext org_mitk_gui_qt_datamanager org_mitk_gui_qt_segmentation org_mitk_gui_qt_volumevisualization org_mitk_gui_qt_diffusionimaging org_mitk_gui_qt_diffusionimaging_connectomics org_mitk_gui_qt_diffusionimaging_fiberfox org_mitk_gui_qt_diffusionimaging_fiberprocessing org_mitk_gui_qt_diffusionimaging_ivim org_mitk_gui_qt_diffusionimaging_odfpeaks org_mitk_gui_qt_diffusionimaging_preprocessing org_mitk_gui_qt_diffusionimaging_reconstruction org_mitk_gui_qt_diffusionimaging_tractography + org_mitk_gui_qt_matchpoint_algorithm_browser + org_mitk_gui_qt_matchpoint_algorithm_control + org_mitk_gui_qt_matchpoint_mapper org_mitk_gui_qt_imagenavigator org_mitk_gui_qt_moviemaker org_mitk_gui_qt_measurementtoolbox org_mitk_gui_qt_basicimageprocessing org_mitk_gui_qt_viewnavigator ) diff --git a/CMake/BuildConfigurations/mitkDiffusion.cmake b/CMake/BuildConfigurations/DiffusionAll.cmake similarity index 94% rename from CMake/BuildConfigurations/mitkDiffusion.cmake rename to CMake/BuildConfigurations/DiffusionAll.cmake index 5e00e6b54b..26f661e2c9 100644 --- a/CMake/BuildConfigurations/mitkDiffusion.cmake +++ b/CMake/BuildConfigurations/DiffusionAll.cmake @@ -1,35 +1,34 @@ -message(STATUS "Configuring MITK Diffusion Release Build") +message(STATUS "Configuring MITK Diffusion with all Plugins") # Enable non-optional external dependencies set(MITK_USE_Vigra ON CACHE BOOL "MITK Use Vigra Library" FORCE) set(MITK_USE_HDF5 ON CACHE BOOL "MITK Use HDF5 Library" FORCE) set(MITK_USE_DCMQI ON CACHE BOOL "" FORCE) set(MITK_USE_MatchPoint ON CACHE BOOL "" FORCE) # Disable all apps but MITK Diffusion set(MITK_BUILD_ALL_APPS OFF CACHE BOOL "Build all MITK applications" FORCE) set(MITK_BUILD_APP_CoreApp OFF CACHE BOOL "Build the MITK CoreApp" FORCE) -set(MITK_BUILD_APP_Workbench OFF CACHE BOOL "Build the MITK Workbench" FORCE) set(MITK_BUILD_APP_Diffusion ON CACHE BOOL "Build MITK Diffusion" FORCE) # Activate Diffusion Mini Apps set(BUILD_DiffusionCoreCmdApps ON CACHE BOOL "Build commandline tools for diffusion" FORCE) set(BUILD_DiffusionFiberProcessingCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber processing" FORCE) set(BUILD_DiffusionFiberfoxCmdApps ON CACHE BOOL "Build commandline tools for diffusion data simulation (Fiberfox)" FORCE) set(BUILD_DiffusionMiscCmdApps ON CACHE BOOL "Build miscellaneous commandline tools for diffusion" FORCE) set(BUILD_DiffusionQuantificationCmdApps ON CACHE BOOL "Build commandline tools for diffusion quantification (IVIM, ADC, ...)" FORCE) set(BUILD_DiffusionTractographyCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography" FORCE) set(BUILD_DiffusionTractographyEvaluationCmdApps ON CACHE BOOL "Build commandline tools for diffusion fiber tractography evaluation" FORCE) set(BUILD_DiffusionConnectomicsCmdApps ON CACHE BOOL "Build commandline tools for diffusion connectomics" FORCE) # Build neither all plugins nor examples set(MITK_BUILD_ALL_PLUGINS OFF CACHE BOOL "Build all MITK plugins" FORCE) set(MITK_BUILD_EXAMPLES OFF CACHE BOOL "Build the MITK examples" FORCE) set(BUILD_TESTING ON CACHE BOOL "Build the MITK tests" FORCE) # Activate in-application help generation set(MITK_DOXYGEN_GENERATE_QCH_FILES ON CACHE BOOL "Use doxygen to generate Qt compressed help files for MITK docs" FORCE) set(BLUEBERRY_USE_QT_HELP ON CACHE BOOL "Enable support for integrating bundle documentation into Qt Help" FORCE) # Disable console window set(MITK_SHOW_CONSOLE_WINDOW ON CACHE BOOL "Use this to enable or disable the console window when starting MITK GUI Applications" FORCE) diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp index d9a719ec3d..5b9795c621 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.cpp @@ -1,437 +1,447 @@ #include "itkFitFibersToImageFilter.h" #include namespace itk{ FitFibersToImageFilter::FitFibersToImageFilter() : m_PeakImage(nullptr) , m_MaskImage(nullptr) , m_FitIndividualFibers(true) , m_GradientTolerance(1e-5) , m_Lambda(0.1) , m_MaxIterations(20) , m_FiberSampling(10) , m_Coverage(0) , m_Overshoot(0) , m_FilterOutliers(true) , m_MeanWeight(1.0) , m_MedianWeight(1.0) , m_MinWeight(1.0) , m_MaxWeight(1.0) , m_Verbose(true) , m_DeepCopy(true) + , m_ResampleFibers(true) , m_NumUnknowns(0) , m_NumResiduals(0) , m_NumCoveredDirections(0) { this->SetNumberOfRequiredOutputs(3); } FitFibersToImageFilter::~FitFibersToImageFilter() { } void FitFibersToImageFilter::GenerateData() { int sz_x = m_PeakImage->GetLargestPossibleRegion().GetSize(0); int sz_y = m_PeakImage->GetLargestPossibleRegion().GetSize(1); int sz_z = m_PeakImage->GetLargestPossibleRegion().GetSize(2); int sz_peaks = m_PeakImage->GetLargestPossibleRegion().GetSize(3)/3 + 1; // +1 for zero - peak int num_voxels = sz_x*sz_y*sz_z; float minSpacing = 1; if(m_PeakImage->GetSpacing()[0]GetSpacing()[1] && m_PeakImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_PeakImage->GetSpacing()[0]; else if (m_PeakImage->GetSpacing()[1] < m_PeakImage->GetSpacing()[2]) minSpacing = m_PeakImage->GetSpacing()[1]; else minSpacing = m_PeakImage->GetSpacing()[2]; m_NumUnknowns = m_Tractograms.size(); if (m_FitIndividualFibers) { m_NumUnknowns = 0; for (unsigned int bundle=0; bundleGetNumFibers(); } - for (unsigned int bundle=0; bundleGetDeepCopy(); - m_Tractograms.at(bundle)->ResampleLinear(minSpacing/m_FiberSampling); - } + if (m_ResampleFibers) + for (unsigned int bundle=0; bundleGetDeepCopy(); + m_Tractograms.at(bundle)->ResampleLinear(minSpacing/m_FiberSampling); + std::cout.rdbuf (old); + } m_NumResiduals = num_voxels * sz_peaks; MITK_INFO << "Num. unknowns: " << m_NumUnknowns; MITK_INFO << "Num. residuals: " << m_NumResiduals; MITK_INFO << "Creating system ..."; vnl_sparse_matrix A; vnl_vector b; A.set_size(m_NumResiduals, m_NumUnknowns); b.set_size(m_NumResiduals); b.fill(0.0); double TD = 0; double FD = 0; m_NumCoveredDirections = 0; unsigned int fiber_count = 0; for (unsigned int bundle=0; bundle polydata = m_Tractograms.at(bundle)->GetFiberPolyData(); for (int i=0; iGetNumFibers(); ++i) { vtkCell* cell = polydata->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (numPoints<2) MITK_INFO << "FIBER WITH ONLY ONE POINT ENCOUNTERED!"; for (int j=0; jGetPoint(j); PointType4 p; p[0]=p1[0]; p[1]=p1[1]; p[2]=p1[2]; p[3]=0; itk::Index<4> idx4; m_PeakImage->TransformPhysicalPointToIndex(p, idx4); itk::Index<3> idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2]; if (!m_PeakImage->GetLargestPossibleRegion().IsInside(idx4) || (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0)) continue; double* p2 = points->GetPoint(j+1); vnl_vector_fixed fiber_dir; fiber_dir[0] = p[0]-p2[0]; fiber_dir[1] = p[1]-p2[1]; fiber_dir[2] = p[2]-p2[2]; fiber_dir.normalize(); double w = 1; int peak_id = sz_peaks-1; vnl_vector_fixed odf_peak = GetClosestPeak(idx4, m_PeakImage, fiber_dir, peak_id, w); float peak_mag = odf_peak.magnitude(); int x = idx4[0]; int y = idx4[1]; int z = idx4[2]; unsigned int linear_index = x + sz_x*y + sz_x*sz_y*z + sz_x*sz_y*sz_z*peak_id; if (b[linear_index] == 0 && peak_id<3) { m_NumCoveredDirections++; FD += peak_mag; } TD += w; if (m_FitIndividualFibers) { b[linear_index] = peak_mag; A.put(linear_index, fiber_count, A.get(linear_index, fiber_count) + w); } else { b[linear_index] = peak_mag; A.put(linear_index, bundle, A.get(linear_index, bundle) + w); } } ++fiber_count; } } TD /= (m_NumCoveredDirections*fiber_count); FD /= m_NumCoveredDirections; A /= TD; b *= 100.0/FD; // times 100 because we want to avoid too small values for computational reasons double init_lambda = fiber_count; // initialization for lambda estimation itk::TimeProbe clock; clock.Start(); VnlCostFunction cost(m_NumUnknowns); cost.SetProblem(A, b, init_lambda); m_Weights.set_size(m_NumUnknowns); // m_Weights.fill( TD/100.0 * FD/2.0 ); m_Weights.fill( 0.0 ); vnl_lbfgsb minimizer(cost); vnl_vector l; l.set_size(m_NumUnknowns); l.fill(0); vnl_vector bound_selection; bound_selection.set_size(m_NumUnknowns); bound_selection.fill(1); minimizer.set_bound_selection(bound_selection); minimizer.set_lower_bound(l); minimizer.set_projected_gradient_tolerance(m_GradientTolerance); MITK_INFO << "Estimating regularization"; minimizer.set_trace(false); minimizer.set_max_function_evals(2); minimizer.minimize(m_Weights); vnl_vector dx; dx.set_size(m_NumUnknowns); dx.fill(0.0); cost.grad_regu_localMSE(m_Weights, dx); double r = dx.magnitude()/m_Weights.magnitude(); cost.m_Lambda *= m_Lambda*55.0/r; MITK_INFO << r << " - " << m_Lambda*55.0/r; if (cost.m_Lambda>10e7) { MITK_INFO << "Regularization estimation failed. Using default value."; cost.m_Lambda = fiber_count; } MITK_INFO << "Using regularization factor of " << cost.m_Lambda << " (λ: " << m_Lambda << ")"; MITK_INFO << "Fitting fibers"; minimizer.set_trace(m_Verbose); minimizer.set_max_function_evals(m_MaxIterations); minimizer.minimize(m_Weights); std::vector< double > weights; if (m_FilterOutliers) { for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); MITK_INFO << "Setting upper weight bound to " << weights.at(m_NumUnknowns*0.99); vnl_vector u; u.set_size(m_NumUnknowns); u.fill(weights.at(m_NumUnknowns*0.99)); minimizer.set_upper_bound(u); bound_selection.fill(2); minimizer.set_bound_selection(bound_selection); minimizer.minimize(m_Weights); weights.clear(); } for (auto w : m_Weights) weights.push_back(w); std::sort(weights.begin(), weights.end()); m_MeanWeight = m_Weights.mean(); m_MedianWeight = weights.at(m_NumUnknowns*0.5); m_MinWeight = weights.at(0); m_MaxWeight = weights.at(m_NumUnknowns-1); MITK_INFO << "*************************"; MITK_INFO << "Weight statistics"; MITK_INFO << "Mean: " << m_MeanWeight; MITK_INFO << "Median: " << m_MedianWeight; MITK_INFO << "75% quantile: " << weights.at(m_NumUnknowns*0.75); MITK_INFO << "95% quantile: " << weights.at(m_NumUnknowns*0.95); MITK_INFO << "99% quantile: " << weights.at(m_NumUnknowns*0.99); MITK_INFO << "Min: " << m_MinWeight; MITK_INFO << "Max: " << m_MaxWeight; MITK_INFO << "*************************"; MITK_INFO << "NumEvals: " << minimizer.get_num_evaluations(); MITK_INFO << "NumIterations: " << minimizer.get_num_iterations(); MITK_INFO << "Residual cost: " << minimizer.get_end_error(); MITK_INFO << "Final RMS: " << cost.S->get_rms_error(m_Weights); clock.Stop(); int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "Optimization took " << h << "h, " << m << "m and " << s << "s"; // transform back for peak image creation A *= FD/100.0; b *= FD/100.0; MITK_INFO << "Weighting fibers"; + std::streambuf *old = cout.rdbuf(); // <-- save + std::stringstream ss; + std::cout.rdbuf (ss.rdbuf()); if (m_FitIndividualFibers) { unsigned int fiber_count = 0; for (unsigned int bundle=0; bundleGetNumFibers(); i++) { m_Tractograms.at(bundle)->SetFiberWeight(i, m_Weights[fiber_count]); ++fiber_count; } m_Tractograms.at(bundle)->Compress(0.1); m_Tractograms.at(bundle)->ColorFibersByFiberWeights(false, true); } } else { for (unsigned int i=0; iSetFiberWeights(m_Weights[i]); m_Tractograms.at(i)->Compress(0.1); m_Tractograms.at(i)->ColorFibersByFiberWeights(false, true); } } + std::cout.rdbuf (old); MITK_INFO << "Generating output images ..."; itk::ImageDuplicator< PeakImgType >::Pointer duplicator = itk::ImageDuplicator< PeakImgType >::New(); duplicator->SetInputImage(m_PeakImage); duplicator->Update(); m_UnderexplainedImage = duplicator->GetOutput(); m_UnderexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_UnderexplainedImage); duplicator->Update(); m_OverexplainedImage = duplicator->GetOutput(); m_OverexplainedImage->FillBuffer(0.0); duplicator->SetInputImage(m_OverexplainedImage); duplicator->Update(); m_ResidualImage = duplicator->GetOutput(); m_ResidualImage->FillBuffer(0.0); duplicator->SetInputImage(m_ResidualImage); duplicator->Update(); m_FittedImage = duplicator->GetOutput(); m_FittedImage->FillBuffer(0.0); vnl_vector fitted_b; fitted_b.set_size(b.size()); cost.S->multiply(m_Weights, fitted_b); for (unsigned int r=0; r idx4; unsigned int linear_index = r; idx4[0] = linear_index % sz_x; linear_index /= sz_x; idx4[1] = linear_index % sz_y; linear_index /= sz_y; idx4[2] = linear_index % sz_z; linear_index /= sz_z; int peak_id = linear_index % sz_peaks; if (peak_id peak_dir; idx4[3] = peak_id*3; peak_dir[0] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[1] = m_PeakImage->GetPixel(idx4); idx4[3] += 1; peak_dir[2] = m_PeakImage->GetPixel(idx4); peak_dir.normalize(); peak_dir *= fitted_b[r]; idx4[3] = peak_id*3; m_FittedImage->SetPixel(idx4, peak_dir[0]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[1]); idx4[3] += 1; m_FittedImage->SetPixel(idx4, peak_dir[2]); } } FD = 0; m_Coverage = 0; m_Overshoot = 0; itk::Index<4> idx4; for (idx4[0]=0; idx4[0] idx3; idx3[0] = idx4[0]; idx3[1] = idx4[1]; idx3[2] = idx4[2]; if (m_MaskImage.IsNotNull() && m_MaskImage->GetPixel(idx3)==0) continue; vnl_vector_fixed peak_dir; vnl_vector_fixed fitted_dir; vnl_vector_fixed overshoot_dir; for (idx4[3]=0; idx4[3]<(itk::IndexValueType)m_PeakImage->GetLargestPossibleRegion().GetSize(3); ++idx4[3]) { peak_dir[idx4[3]%3] = m_PeakImage->GetPixel(idx4); fitted_dir[idx4[3]%3] = m_FittedImage->GetPixel(idx4); m_ResidualImage->SetPixel(idx4, m_PeakImage->GetPixel(idx4) - m_FittedImage->GetPixel(idx4)); if (idx4[3]%3==2) { FD += peak_dir.magnitude(); itk::Index<4> tidx= idx4; if (peak_dir.magnitude()>fitted_dir.magnitude()) { m_Coverage += fitted_dir.magnitude(); m_UnderexplainedImage->SetPixel(tidx, peak_dir[2]-fitted_dir[2]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[1]-fitted_dir[1]); tidx[3]--; m_UnderexplainedImage->SetPixel(tidx, peak_dir[0]-fitted_dir[0]); } else { overshoot_dir[0] = fitted_dir[0]-peak_dir[0]; overshoot_dir[1] = fitted_dir[1]-peak_dir[1]; overshoot_dir[2] = fitted_dir[2]-peak_dir[2]; m_Coverage += peak_dir.magnitude(); m_Overshoot += overshoot_dir.magnitude(); m_OverexplainedImage->SetPixel(tidx, overshoot_dir[2]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[1]); tidx[3]--; m_OverexplainedImage->SetPixel(tidx, overshoot_dir[0]); } } } } m_Coverage = m_Coverage/FD; m_Overshoot = m_Overshoot/FD; - MITK_INFO << std::fixed << "Coverage: " << setprecision(1) << 100.0*m_Coverage << "%"; - MITK_INFO << std::fixed << "Overshoot: " << setprecision(1) << 100.0*m_Overshoot << "%"; + MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*m_Coverage << "%"; + MITK_INFO << std::fixed << "Overshoot: " << setprecision(2) << 100.0*m_Overshoot << "%"; } vnl_vector_fixed FitFibersToImageFilter::GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer peak_image , vnl_vector_fixed fiber_dir, int& id, double& w ) { int m_NumDirs = peak_image->GetLargestPossibleRegion().GetSize()[3]/3; vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0.8; for (int i=0; i dir; idx[3] = i*3; dir[0] = peak_image->GetPixel(idx); idx[3] += 1; dir[1] = peak_image->GetPixel(idx); idx[3] += 1; dir[2] = peak_image->GetPixel(idx); float mag = dir.magnitude(); if (magangle) { angle = fabs(a); w = angle; if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= mag; id = i; } } return out_dir; } std::vector FitFibersToImageFilter::GetTractograms() const { return m_Tractograms; } void FitFibersToImageFilter::SetTractograms(const std::vector &tractograms) { m_Tractograms = tractograms; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h index 119d39e57d..16c060af2b 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFitFibersToImageFilter.h @@ -1,262 +1,265 @@ #ifndef __itkFitFibersToImageFilter_h__ #define __itkFitFibersToImageFilter_h__ // MITK #include #include #include #include #include #include #include #include #include namespace itk{ /** * \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). */ class FitFibersToImageFilter : public ImageSource< mitk::PeakImage::ItkPeakImageType > { public: typedef FitFibersToImageFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; typedef itk::Point PointType4; typedef mitk::PeakImage::ItkPeakImageType PeakImgType; typedef itk::Image UcharImgType; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( FitFibersToImageFilter, ImageSource ) itkSetMacro( PeakImage, PeakImgType::Pointer) itkGetMacro( PeakImage, PeakImgType::Pointer) itkSetMacro( MaskImage, UcharImgType::Pointer) itkGetMacro( MaskImage, UcharImgType::Pointer) itkSetMacro( FitIndividualFibers, bool) itkGetMacro( FitIndividualFibers, bool) itkSetMacro( GradientTolerance, double) itkGetMacro( GradientTolerance, double) itkSetMacro( Lambda, double) itkGetMacro( Lambda, double) itkSetMacro( MaxIterations, int) itkGetMacro( MaxIterations, int) itkSetMacro( FiberSampling, float) itkGetMacro( FiberSampling, float) itkSetMacro( FilterOutliers, bool) itkGetMacro( FilterOutliers, bool) itkSetMacro( Verbose, bool) itkGetMacro( Verbose, bool) itkSetMacro( DeepCopy, bool) itkGetMacro( DeepCopy, bool) + itkSetMacro( ResampleFibers, bool) + itkGetMacro( ResampleFibers, bool) itkGetMacro( Weights, vnl_vector) itkGetMacro( FittedImage, PeakImgType::Pointer) itkGetMacro( ResidualImage, PeakImgType::Pointer) itkGetMacro( OverexplainedImage, PeakImgType::Pointer) itkGetMacro( UnderexplainedImage, PeakImgType::Pointer) itkGetMacro( Coverage, double) itkGetMacro( Overshoot, double) itkGetMacro( MeanWeight, double) itkGetMacro( MedianWeight, double) itkGetMacro( MinWeight, double) itkGetMacro( MaxWeight, double) itkGetMacro( NumUnknowns, unsigned int) itkGetMacro( NumResiduals, unsigned int) itkGetMacro( NumCoveredDirections, unsigned int) void SetTractograms(const std::vector &tractograms); void GenerateData() override; std::vector GetTractograms() const; protected: FitFibersToImageFilter(); virtual ~FitFibersToImageFilter(); vnl_vector_fixed GetClosestPeak(itk::Index<4> idx, PeakImgType::Pointer m_PeakImage , vnl_vector_fixed fiber_dir, int& id, double& w ); std::vector< mitk::FiberBundle::Pointer > m_Tractograms; PeakImgType::Pointer m_PeakImage; UcharImgType::Pointer m_MaskImage; bool m_FitIndividualFibers; double m_GradientTolerance; double m_Lambda; int m_MaxIterations; float m_FiberSampling; double m_Coverage; double m_Overshoot; bool m_FilterOutliers; double m_MeanWeight; double m_MedianWeight; double m_MinWeight; double m_MaxWeight; bool m_Verbose; bool m_DeepCopy; + bool m_ResampleFibers; unsigned int m_NumUnknowns; unsigned int m_NumResiduals; unsigned int m_NumCoveredDirections; // output vnl_vector m_Weights; PeakImgType::Pointer m_UnderexplainedImage; PeakImgType::Pointer m_OverexplainedImage; PeakImgType::Pointer m_ResidualImage; PeakImgType::Pointer m_FittedImage; }; } class VnlCostFunction : public vnl_cost_function { public: vnl_sparse_matrix_linear_system< double >* S; vnl_sparse_matrix< double > m_A; vnl_sparse_matrix< double > m_A_Ones; // matrix indicating active weights with 1 vnl_vector< double > m_b; double m_Lambda; // regularization factor vnl_vector row_sums; // number of active weights per row vnl_vector local_weight_means; // mean weight of each row void SetProblem(vnl_sparse_matrix< double >& A, vnl_vector& b, double lambda) { S = new vnl_sparse_matrix_linear_system(A, b); m_A = A; m_b = b; m_Lambda = lambda; m_A_Ones.set_size(m_A.rows(), m_A.cols()); m_A.reset(); while (m_A.next()) m_A_Ones.put(m_A.getrow(), m_A.getcolumn(), 1); unsigned int N = m_b.size(); vnl_vector ones; ones.set_size(dim); ones.fill(1.0); row_sums.set_size(N); m_A_Ones.mult(ones, row_sums); local_weight_means.set_size(N); } VnlCostFunction(const int NumVars) : vnl_cost_function(NumVars) { } void regu_MSE(vnl_vector const &x, double& cost) { double mean = x.mean(); vnl_vector tx = x-mean; cost += m_Lambda*1e8*tx.squared_magnitude()/x.size(); } void regu_MSM(vnl_vector const &x, double& cost) { cost += m_Lambda*1e8*x.squared_magnitude()/x.size(); } void regu_localMSE(vnl_vector const &x, double& cost) { m_A_Ones.mult(x, local_weight_means); local_weight_means = element_quotient(local_weight_means, row_sums); m_A_Ones.reset(); double regu = 0; while (m_A_Ones.next()) { double d = 0; if (x[m_A_Ones.getcolumn()]>local_weight_means[m_A_Ones.getrow()]) d = std::exp(x[m_A_Ones.getcolumn()]) - std::exp(local_weight_means[m_A_Ones.getrow()]); else d = x[m_A_Ones.getcolumn()] - local_weight_means[m_A_Ones.getrow()]; regu += d*d; } cost += m_Lambda*regu/dim; } void grad_regu_MSE(vnl_vector const &x, vnl_vector &dx) { double mean = x.mean(); vnl_vector tx = x-mean; vnl_vector tx2(dim, 0.0); vnl_vector h(dim, 1.0); for (int c=0; c const &x, vnl_vector &dx) { dx += m_Lambda*1e8*2.0*x/dim; } void grad_regu_localMSE(vnl_vector const &x, vnl_vector &dx) { m_A_Ones.mult(x, local_weight_means); local_weight_means = element_quotient(local_weight_means, row_sums); vnl_vector exp_x = x.apply(std::exp); vnl_vector exp_means = local_weight_means.apply(std::exp); vnl_vector tdx(dim, 0); m_A_Ones.reset(); while (m_A_Ones.next()) { int c = m_A_Ones.getcolumn(); int r = m_A_Ones.getrow(); if (x[c]>local_weight_means[r]) tdx[c] += exp_x[c] * ( exp_x[c] - exp_means[r] ); else tdx[c] += x[c] - local_weight_means[r]; } dx += tdx*2.0*m_Lambda/dim; } double f(vnl_vector const &x) { double cost = S->get_rms_error(x); cost *= cost; regu_localMSE(x, cost); return cost; } void gradf(vnl_vector const &x, vnl_vector &dx) { dx.fill(0.0); unsigned int N = m_b.size(); vnl_vector d; d.set_size(N); S->multiply(x,d); d -= m_b; S->transpose_multiply(d, dx); dx *= 2.0/N; grad_regu_localMSE(x,dx); } }; #ifndef ITK_MANUAL_INSTANTIATION #include "itkFitFibersToImageFilter.cpp" #endif #endif // __itkFitFibersToImageFilter_h__ diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp index d33bf2cf32..8092e94a0d 100755 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp @@ -1,2335 +1,2383 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkFiberBundle.h" #include #include #include #include "mitkImagePixelReadAccessor.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs"; using namespace std; mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData ) : m_NumFibers(0) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != nullptr) m_FiberPolyData = fiberPolyData; this->UpdateFiberGeometry(); this->GenerateFiberIds(); this->ColorFibersByOrientation(); } mitk::FiberBundle::~FiberBundle() { } mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy() { mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); newFib->SetFiberColors(this->m_FiberColors); newFib->SetFiberWeights(this->m_FiberWeights); return newFib; } vtkSmartPointer mitk::FiberBundle::GeneratePolyDataByIds(std::vector fiberIds) { vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); auto finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*finIt < 0 || *finIt>GetNumFibers()){ MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; break; } vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); vtkSmartPointer fibPoints = fiber->GetPoints(); vtkSmartPointer newFiber = vtkSmartPointer::New(); newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); for(int i=0; iGetNumberOfPoints(); i++) { newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); } newLineSet->InsertNextCell(newFiber); ++finIt; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib) { if (fib==nullptr) { MITK_WARN << "trying to call AddBundle with nullptr argument"; return nullptr; } MITK_INFO << "Adding fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); unsigned int counter = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // add new fiber bundle for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, fib->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // subtract two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib) { MITK_INFO << "Subtracting fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector< std::vector< itk::Point > > points1; for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = GetItkPoint(points->GetPoint(0)); itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); points1.push_back( {start, end} ); } std::vector< std::vector< itk::Point > > points2; for( int i=0; iGetNumFibers(); i++ ) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = GetItkPoint(points->GetPoint(0)); itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); points2.push_back( {start, end} ); } int progress = 0; std::vector< int > ids; #pragma omp parallel for for (int i=0; i<(int)points1.size(); i++) { #pragma omp critical { progress++; std::cout << (int)(100*(float)progress/points1.size()) << "%" << '\r'; cout.flush(); } bool match = false; for (unsigned int j=0; jGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(points->GetPoint(j)); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if(vNewLines->GetNumberOfCells()==0) return nullptr; // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle return mitk::FiberBundle::New(vNewPolyData); } itk::Point mitk::FiberBundle::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } /* * set PolyData (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { if (fiberPD == nullptr) this->m_FiberPolyData = vtkSmartPointer::New(); else m_FiberPolyData->DeepCopy(fiberPD); m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (updateGeometry) UpdateFiberGeometry(); GenerateFiberIds(); ColorFibersByOrientation(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundle::GetFiberPolyData() const { return m_FiberPolyData; } void mitk::FiberBundle::ColorFibersByOrientation() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also PolyData needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= vtkPoints* extrPoints = nullptr; extrPoints = m_FiberPolyData->GetPoints(); int numOfPoints = 0; if (extrPoints!=nullptr) numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = 4; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(numOfPoints * componentSize); m_FiberColors->SetNumberOfComponents(componentSize); m_FiberColors->SetName("FIBER_COLORS"); int numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) return; /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); /* single fiber checkpoints: is number of points valid */ if (pointsPerFiber > 1) { /* operate on points of single fiber */ for (int i=0; i 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); rgba[3] = (unsigned char) (255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); rgba[3] = (unsigned char) (255.0); } m_FiberColors->InsertTypedTuple(idList[i], rgba); } } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; } else { MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByCurvature(bool, bool normalize) { double window = 5; //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = 4; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize); m_FiberColors->SetNumberOfComponents(componentSize); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); vector< double > values; double min = 1; double max = 0; MITK_INFO << "Coloring fibers by curvature"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures for (int j=0; j > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0); while(dist1) { double p1[3]; points->GetPoint(c-1, p1); double p2[3]; points->GetPoint(c, p2); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==j) meanV += v; c--; } c = j; dist = 0; while(distGetPoint(c, p1); double p2[3]; points->GetPoint(c+1, p2); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==j) meanV += v; c++; } meanV.normalize(); double dev = 0; for (unsigned int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); dev = 1.0-dev/180.0; values.push_back(dev); if (devmax) max = dev; } } unsigned int count = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); for (int j=0; j1) dev = 1; lookupTable->GetColor(dev, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTypedTuple(cell->GetPointId(j), rgba); count++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray) { for(long i=0; iGetNumberOfTuples(); i++) { double faValue = FAValArray->GetValue(i); faValue = faValue * 255.0; m_FiberColors->SetComponent(i,3, (unsigned char) faValue ); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ResetFiberOpacity() { for(long i=0; iGetNumberOfTuples(); i++) m_FiberColors->SetComponent(i,3, 255.0 ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize) { mitkPixelTypeMultiplex3( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } template void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); unsigned char rgba[4] = {0,0,0,0}; vtkPoints* pointSet = m_FiberPolyData->GetPoints(); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); double min = 9999999; double max = -9999999; for(long i=0; iGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; double pixelValue = readimage.GetPixelByWorldCoordinates(px); if (pixelValue>max) max = pixelValue; if (pixelValueGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; double pixelValue = readimage.GetPixelByWorldCoordinates(px); if (normalize) pixelValue = (pixelValue-min)/(max-min); else if (pixelValue>1) pixelValue = 1; double color[3]; lookupTable->GetColor(1-pixelValue, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); if (opacity) rgba[3] = (unsigned char) (255.0 * pixelValue); else rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByFiberWeights(bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); unsigned char rgba[4] = {0,0,0,0}; unsigned int counter = 0; float max = -999999; float min = 999999; for (int i=0; iGetFiberWeight(i); if (weight>max) max = weight; if (weightGetCell(i); int numPoints = cell->GetNumberOfPoints(); double weight = this->GetFiberWeight(i); for (int j=0; j1) v = 1; double color[3]; lookupTable->GetColor(1-v, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); if (opacity) rgba[3] = (unsigned char) (255.0 * v); else rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTypedTuple(counter, rgba); counter++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); unsigned char rgba[4] = {0,0,0,0}; for(long i=0; iGetNumberOfPoints(); ++i) { rgba[0] = (unsigned char) r; rgba[1] = (unsigned char) g; rgba[2] = (unsigned char) b; rgba[3] = (unsigned char) alpha; m_FiberColors->InsertTypedTuple(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::GenerateFiberIds() { if (m_FiberPolyData == nullptr) return; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInputData(m_FiberPolyData); idFiberFilter->CellIdsOn(); // idFiberFilter->PointIdsOn(); // point id's are not needed idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); idFiberFilter->FieldDataOn(); idFiberFilter->Update(); m_FiberIdDataSet = idFiberFilter->GetOutput(); } +float mitk::FiberBundle::GetOverlap(ItkUcharImgType* mask, bool do_resampling) +{ + vtkSmartPointer PolyData = m_FiberPolyData; + mitk::FiberBundle::Pointer fibCopy = this; + if (do_resampling) + { + float minSpacing = 1; + if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) + minSpacing = mask->GetSpacing()[0]; + else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) + minSpacing = mask->GetSpacing()[1]; + else + minSpacing = mask->GetSpacing()[2]; + + fibCopy = this->GetDeepCopy(); + fibCopy->ResampleLinear(minSpacing/5); + PolyData = fibCopy->GetFiberPolyData(); + } + + MITK_INFO << "Calculating overlap"; + int inside = 0; + int outside = 0; + boost::progress_display disp(m_NumFibers); + for (int i=0; iGetCell(i); + int numPoints = cell->GetNumberOfPoints(); + vtkPoints* points = cell->GetPoints(); + + for (int j=0; jGetPoint(j); + itk::Point itkP; + itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; + itk::Index<3> idx; + mask->TransformPhysicalPointToIndex(itkP, idx); + + if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 ) + inside++; + else + outside++; + } + } + + return (float)inside/(inside+outside); +} + mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds, float fraction, bool do_resampling) { vtkSmartPointer PolyData = m_FiberPolyData; mitk::FiberBundle::Pointer fibCopy = this; if (anyPoint && do_resampling) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; fibCopy = this->GetDeepCopy(); fibCopy->ResampleLinear(minSpacing/5); PolyData = fibCopy->GetFiberPolyData(); } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); std::vector< float > new_weights; MITK_INFO << "Extracting fibers with mask image"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cellOriginal = m_FiberPolyData->GetCell(i); int numPointsOriginal = cellOriginal->GetNumberOfPoints(); vtkPoints* pointsOriginal = cellOriginal->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); if (numPoints>1 && numPointsOriginal) { if (anyPoint) { int inside = 0; int outside = 0; if (!invert) { for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 ) { inside++; if (fraction==0) break; } else outside++; } float current_fraction = 0.0; if (inside+outside>0) current_fraction = (float)inside/(inside+outside); if (current_fraction>fraction) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else { bool includeFiber = true; for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) ) { inside++; includeFiber = false; break; } else outside++; } if (includeFiber) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } else { double* start = pointsOriginal->GetPoint(0); itk::Point itkStart; itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; itk::Index<3> idxStart; mask->TransformPhysicalPointToIndex(itkStart, idxStart); double* end = pointsOriginal->GetPoint(numPointsOriginal-1); itk::Point itkEnd; itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; itk::Index<3> idxEnd; mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); if (invert) { if (bothEnds) { if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else { if (bothEnds) { if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } } if (container->GetNumberOfPoints()>0) { new_weights.push_back(fibCopy->GetFiberWeight(i)); vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer newfib = mitk::FiberBundle::New(newPolyData); for (unsigned int i=0; iSetFiberWeight(i, new_weights.at(i)); return newfib; } mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); fibCopy->ResampleLinear(minSpacing/10); vtkSmartPointer PolyData =fibCopy->GetFiberPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Cutting fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); if (numPoints>1) { int newNumPoints = 0; for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if (newNumPoints>0) { vtkNewCells->InsertNextCell(container); newNumPoints = 0; container = vtkSmartPointer::New(); } } if (newNumPoints>0) vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData); newFib->Compress(0.1); return newFib; } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage) { if (roi==nullptr || !(dynamic_cast(roi->GetData()) || dynamic_cast(roi->GetData())) ) return nullptr; std::vector tmp = ExtractFiberIdSubset(roi, storage); if (tmp.size()<=0) return mitk::FiberBundle::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); return mitk::FiberBundle::New(pTmp); } std::vector mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage) { std::vector result; if (roi==nullptr || roi->GetData()==nullptr) return result; mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); if (!pfc.IsNull()) // handle composite { DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); if (children->size()==0) return result; switch (pfc->getOperationType()) { case 0: // AND { MITK_INFO << "AND"; result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(std::min(result.size(),inRoi.size())); it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( it - rest.begin() ); result = rest; } break; } case 1: // OR { MITK_INFO << "OR"; result = ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { it = result.end(); std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); result.insert(it, inRoi.begin(), inRoi.end()); } // remove duplicates sort(result.begin(), result.end()); it = unique(result.begin(), result.end()); result.resize( it - result.begin() ); break; } case 2: // NOT { MITK_INFO << "NOT"; for(long i=0; iGetNumFibers(); i++) result.push_back(i); std::vector::iterator it; for (unsigned int i=0; iSize(); ++i) { std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(result.size()-inRoi.size()); it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( it - rest.begin() ); result = rest; } break; } } } else if ( dynamic_cast(roi->GetData()) ) // actual extraction { if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarPoly = dynamic_cast(roi->GetData()); //create vtkPolygon using controlpoints from planarFigure polygon vtkSmartPointer polygonVtk = vtkSmartPointer::New(); for (unsigned int i=0; iGetNumberOfControlPoints(); ++i) { itk::Point p = planarPoly->GetWorldControlPoint(i); vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] ); polygonVtk->GetPointIds()->InsertNextId(id); } MITK_INFO << "Extracting with polygon"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); double tolerance = 0.001; // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection double pcoords[3] = {0,0,0}; int subId = 0; int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId); if (iD!=0) { result.push_back(i); break; } } } } else if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(roi->GetData()); Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal(); planeNormal.Normalize(); //calculate circle radius mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint double radius = V1w.EuclideanDistanceTo(V2w); radius *= radius; MITK_INFO << "Extracting with circle"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x); if (iD!=0) { double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]); if( dist <= radius) { result.push_back(i); break; } } } } } return result; } return result; } void mitk::FiberBundle::UpdateFiberGeometry() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInputData(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = m_FiberPolyData->GetNumberOfCells(); if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints()) this->ColorFibersByOrientation(); if (m_FiberWeights->GetNumberOfValues()!=m_NumFibers) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberWeights->SetNumberOfValues(m_NumFibers); this->SetFiberWeights(1); } if (m_NumFibers<=0) // no fibers present; apply default geometry { m_MinFiberLength = 0; m_MaxFiberLength = 0; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(false); float b[] = {0, 1, 0, 1, 0, 1}; geometry->SetFloatBounds(b); SetGeometry(geometry); return; } double b[6]; m_FiberPolyData->GetBounds(b); // calculate statistics for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); length += dist; } m_FiberLengths.push_back(length); m_MeanFiberLength += length; if (i==0) { m_MinFiberLength = length; m_MaxFiberLength = length; } else { if (lengthm_MaxFiberLength) m_MaxFiberLength = length; } } m_MeanFiberLength /= m_NumFibers; std::vector< float > sortedLengths = m_FiberLengths; std::sort(sortedLengths.begin(), sortedLengths.end()); for (int i=0; i1) m_LengthStDev /= (m_NumFibers-1); else m_LengthStDev = 0; m_LengthStDev = std::sqrt(m_LengthStDev); m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetFloatBounds(b); this->SetGeometry(geometry); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const { return m_FiberWeights->GetValue(fiber); } void mitk::FiberBundle::SetFiberWeights(float newWeight) { for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, newWeight); } void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer weights) { if (m_NumFibers!=weights->GetNumberOfValues()) { MITK_INFO << "Weights array not equal to number of fibers! " << weights->GetNumberOfValues() << " vs " << m_NumFibers; return; } for (int i=0; iGetNumberOfValues(); i++) m_FiberWeights->SetValue(i, weights->GetValue(i)); m_FiberWeights->SetName("FIBER_WEIGHTS"); } void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight) { m_FiberWeights->SetValue(fiber, weight); } void mitk::FiberBundle::SetFiberColors(vtkSmartPointer fiberColors) { for(long i=0; iGetNumberOfPoints(); ++i) { unsigned char source[4] = {0,0,0,0}; fiberColors->GetTypedTuple(i, source); unsigned char target[4] = {0,0,0,0}; target[0] = source[0]; target[1] = source[1]; target[2] = source[2]; target[3] = source[3]; m_FiberColors->InsertTypedTuple(i, target); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; m = rot*m; return m; } itk::Point mitk::FiberBundle::TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); point[0] -= center[0]; point[1] -= center[1]; point[2] -= center[2]; point = rot*point; point[0] += center[0]+tx; point[1] += center[1]+ty; point[2] += center[2]+tz; itk::Point out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2]; return out; } void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rot*dir; dir[0] += center[0]+tx; dir[1] += center[1]+ty; dir[2] += center[2]+tz; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z) { x = x*M_PI/180; y = y*M_PI/180; z = z*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rotZ*rotY*rotX*dir; dir[0] += center[0]; dir[1] += center[1]; dir[2] += center[2]; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter) { MITK_INFO << "Scaling fibers"; boost::progress_display disp(m_NumFibers); mitk::BaseGeometry* geom = this->GetGeometry(); mitk::Point3D c = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); if (subtractCenter) { p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; } p[0] *= x; p[1] *= y; p[2] *= z; if (subtractCenter) { p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; } vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TranslateFibers(double x, double y, double z) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[0] += x; p[1] += y; p[2] += z; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::MirrorFibers(unsigned int axis) { if (axis>2) return; MITK_INFO << "Mirroring fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[axis] = -p[axis]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RemoveDir(vnl_vector_fixed dir, double threshold) { dir.normalize(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); bool discard = false; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); vnl_vector_fixed< double, 3 > v1; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; if (v1.magnitude()>0.001) { v1.normalize(); if (fabs(dot_product(v1,dir))>threshold) { discard = true; break; } } } if (!discard) { for (int j=0; jGetPoint(j, p1); vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); // UpdateColorCoding(); // UpdateFiberGeometry(); } bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double p3[3]; points->GetPoint(j+2, p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; v2[0] = p3[0]-p2[0]; v2[1] = p3[1]-p2[1]; v2[2] = p3[2]-p2[2]; v3[0] = p1[0]-p3[0]; v3[1] = p1[1]-p3[1]; v3[2] = p1[2]-p3[2]; float a = v1.magnitude(); float b = v2.magnitude(); float c = v3.magnitude(); float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); if (deleteFibers && rInsertNextCell(container); container = vtkSmartPointer::New(); } else if (j==numPoints-3) { id = vtkNewPoints->InsertNextPoint(p2); container->GetPointIds()->InsertNextId(id); id = vtkNewPoints->InsertNextPoint(p3); container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); } } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM) { MITK_INFO << "Removing short fibers"; if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers { MITK_WARN << "Process aborted. No fibers would be left!"; return false; } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); float min = m_MaxFiberLength; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); if (m_FiberLengths.at(i)GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Removing long fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias ) { if (pointDistance<=0) return; vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::New(); //in smoothpoints the interpolated points representing a fiber are stored. //in vtkcells all polylines are stored, actually all id's of them are stored vtkSmartPointer vtkSmoothCells = vtkSmartPointer::New(); //cellcontainer for smoothed lines vtkIdType pointHelperCnt = 0; MITK_INFO << "Smoothing fibers"; vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); boost::progress_display disp(m_NumFibers); #pragma omp parallel for for (int i=0; i newPoints = vtkSmartPointer::New(); float length = 0; float weight = 1; #pragma omp critical { length = m_FiberLengths.at(i); weight = m_FiberWeights->GetValue(i); ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jInsertNextPoint(points->GetPoint(j)); } int sampling = std::ceil(length/pointDistance); vtkSmartPointer xSpline = vtkSmartPointer::New(); vtkSmartPointer ySpline = vtkSmartPointer::New(); vtkSmartPointer zSpline = vtkSmartPointer::New(); xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); vtkSmartPointer spline = vtkSmartPointer::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(newPoints); vtkSmartPointer functionSource = vtkSmartPointer::New(); functionSource->SetParametricFunction(spline); functionSource->SetUResolution(sampling); functionSource->SetVResolution(sampling); functionSource->SetWResolution(sampling); functionSource->Update(); vtkPolyData* outputFunction = functionSource->GetOutput(); vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber vtkSmartPointer smoothLine = vtkSmartPointer::New(); smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints()); #pragma omp critical { for (int j=0; jGetNumberOfPoints(); j++) { smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt); vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); } newFiberWeights->SetValue(vtkSmoothCells->GetNumberOfCells(), weight); vtkSmoothCells->InsertNextCell(smoothLine); pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints(); } } SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ResampleSpline(float pointDistance) { ResampleSpline(pointDistance, 0, 0, 0 ); } unsigned long mitk::FiberBundle::GetNumberOfPoints() const { unsigned long points = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); points += cell->GetNumberOfPoints(); } return points; } void mitk::FiberBundle::Compress(float error) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Compressing fibers"; unsigned long numRemovedPoints = 0; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } // calculate curvatures int numPoints = vertices.size(); std::vector< int > removedPoints; removedPoints.resize(numPoints, 0); removedPoints[0]=-1; removedPoints[numPoints-1]=-1; vtkSmartPointer container = vtkSmartPointer::New(); int remCounter = 0; bool pointFound = true; while (pointFound) { pointFound = false; double minError = error; int removeIndex = -1; for (unsigned int j=0; j candV = vertices.at(j); int validP = -1; vnl_vector_fixed< double, 3 > pred; for (int k=j-1; k>=0; k--) if (removedPoints[k]<=0) { pred = vertices.at(k); validP = k; break; } int validS = -1; vnl_vector_fixed< double, 3 > succ; for (int k=j+1; k=0 && validS>=0) { double a = (candV-pred).magnitude(); double b = (candV-succ).magnitude(); double c = (pred-succ).magnitude(); double s=0.5*(a+b+c); double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c))); if (hcInsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); numRemovedPoints += remCounter; vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { MITK_INFO << "Removed points: " << numRemovedPoints; SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } void mitk::FiberBundle::ResampleLinear(double pointDistance) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Resampling fibers (linear)"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); container->GetPointIds()->InsertNextId(id); } for (unsigned int j=1; j vec = vertices.at(j) - lastV; double new_dist = vec.magnitude(); if (new_dist >= pointDistance) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-pointDistance <= mitk::eps ) { vec.normalize(); newV += vec * pointDistance; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if (j==vertices.size()-1 && new_dist>0.0001) { #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } // reapply selected colorcoding in case PolyData structure has changed bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps) { if (fib==nullptr) { MITK_INFO << "Reference bundle is nullptr!"; return false; } if (m_NumFibers!=fib->GetNumFibers()) { MITK_INFO << "Unequal number of fibers!"; MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); return false; } for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i); int numPoints2 = cell2->GetNumberOfPoints(); vtkPoints* points2 = cell2->GetPoints(); if (numPoints2!=numPoints) { MITK_INFO << "Unequal number of points in fiber " << i << "!"; MITK_INFO << numPoints2 << " vs. " << numPoints; return false; } for (int j=0; jGetPoint(j); double* p2 = points2->GetPoint(j); if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps) { MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!"; MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2]; MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2]; return false; } } } return true; } void mitk::FiberBundle::PrintSelf(std::ostream &os, itk::Indent indent) const { os << indent << this->GetNameOfClass() << ":\n"; os << indent << "Number of fibers: " << this->GetNumFibers() << std::endl; os << indent << "Min. fiber length: " << this->GetMinFiberLength() << std::endl; os << indent << "Max. fiber length: " << this->GetMaxFiberLength() << std::endl; os << indent << "Mean fiber length: " << this->GetMeanFiberLength() << std::endl; os << indent << "Median fiber length: " << this->GetMedianFiberLength() << std::endl; os << indent << "STDEV fiber length: " << this->GetLengthStDev() << std::endl; os << indent << "Number of points: " << this->GetNumberOfPoints() << std::endl; Superclass::PrintSelf(os, indent); } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundle::UpdateOutputInformation() { } void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundle::VerifyRequestedRegion() { return true; } void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* ) { } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h index a49a44e2dd..862cfd00a3 100644 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.h @@ -1,179 +1,176 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _MITK_FiberBundle_H #define _MITK_FiberBundle_H //includes for MITK datastructure #include #include #include #include #include #include #include //includes storing fiberdata #include #include #include #include #include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class MITKFIBERTRACKING_EXPORT FiberBundle : public BaseData { public: typedef itk::Image ItkUcharImgType; // fiber colorcodings static const char* FIBER_ID_ARRAY; virtual void UpdateOutputInformation() override; virtual void SetRequestedRegionToLargestPossibleRegion() override; virtual bool RequestedRegionIsOutsideOfTheBufferedRegion() override; virtual bool VerifyRequestedRegion() override; virtual void SetRequestedRegion(const itk::DataObject*) override; mitkClassMacro( FiberBundle, BaseData ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) mitkNewMacro1Param(Self, vtkSmartPointer) // custom constructor // colorcoding related methods void ColorFibersByFiberWeights(bool opacity, bool normalize); void ColorFibersByCurvature(bool opacity, bool normalize); void ColorFibersByScalarMap(mitk::Image::Pointer, bool opacity, bool normalize); template void ColorFibersByScalarMap(const mitk::PixelType pixelType, mitk::Image::Pointer, bool opacity, bool normalize); void ColorFibersByOrientation(); void SetFiberOpacity(vtkDoubleArray *FAValArray); void ResetFiberOpacity(); void SetFiberColors(vtkSmartPointer fiberColors); void SetFiberColors(float r, float g, float b, float alpha=255); vtkSmartPointer GetFiberColors() const { return m_FiberColors; } // fiber compression void Compress(float error = 0.0); // fiber resampling void ResampleSpline(float pointDistance=1); void ResampleSpline(float pointDistance, double tension, double continuity, double bias ); void ResampleLinear(double pointDistance=1); bool RemoveShortFibers(float lengthInMM); bool RemoveLongFibers(float lengthInMM); bool ApplyCurvatureThreshold(float minRadius, bool deleteFibers); void MirrorFibers(unsigned int axis); void RotateAroundAxis(double x, double y, double z); void TranslateFibers(double x, double y, double z); void ScaleFibers(double x, double y, double z, bool subtractCenter=true); void TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz); void RemoveDir(vnl_vector_fixed dir, double threshold); itk::Point TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz); itk::Matrix< double, 3, 3 > TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz); // add/subtract fibers FiberBundle::Pointer AddBundle(FiberBundle* fib); FiberBundle::Pointer SubtractBundle(FiberBundle* fib); // fiber subset extraction FiberBundle::Pointer ExtractFiberSubset(DataNode *roi, DataStorage* storage); std::vector ExtractFiberIdSubset(DataNode* roi, DataStorage* storage); FiberBundle::Pointer ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert=false, bool bothEnds=true, float fraction=0.0, bool do_resampling=true); FiberBundle::Pointer RemoveFibersOutside(ItkUcharImgType* mask, bool invert=false); - - vtkSmartPointer GeneratePolyDataByIds( std::vector ); // TODO: make protected - void GenerateFiberIds(); // TODO: make protected + float GetOverlap(ItkUcharImgType* mask, bool do_resampling); // get/set data vtkSmartPointer GetFiberWeights() const { return m_FiberWeights; } float GetFiberWeight(unsigned int fiber) const; void SetFiberWeights(float newWeight); void SetFiberWeight(unsigned int fiber, float weight); void SetFiberWeights(vtkSmartPointer weights); void SetFiberPolyData(vtkSmartPointer, bool updateGeometry = true); vtkSmartPointer GetFiberPolyData() const; itkGetConstMacro( NumFibers, int) //itkGetMacro( FiberSampling, int) itkGetConstMacro( MinFiberLength, float ) itkGetConstMacro( MaxFiberLength, float ) itkGetConstMacro( MeanFiberLength, float ) itkGetConstMacro( MedianFiberLength, float ) itkGetConstMacro( LengthStDev, float ) itkGetConstMacro( UpdateTime2D, itk::TimeStamp ) itkGetConstMacro( UpdateTime3D, itk::TimeStamp ) void RequestUpdate2D(){ m_UpdateTime2D.Modified(); } void RequestUpdate3D(){ m_UpdateTime3D.Modified(); } void RequestUpdate(){ m_UpdateTime2D.Modified(); m_UpdateTime3D.Modified(); } unsigned long GetNumberOfPoints() const; // copy fiber bundle mitk::FiberBundle::Pointer GetDeepCopy(); // compare fiber bundles bool Equals(FiberBundle* fib, double eps=0.01); itkSetMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) itkGetConstMacro( ReferenceGeometry, mitk::BaseGeometry::Pointer ) protected: FiberBundle( vtkPolyData* fiberPolyData = nullptr ); virtual ~FiberBundle(); - itk::Point GetItkPoint(double point[3]); - - // calculate geometry from fiber extent - void UpdateFiberGeometry(); - - virtual void PrintSelf(std::ostream &os, itk::Indent indent) const override; + vtkSmartPointer GeneratePolyDataByIds( std::vector ); + void GenerateFiberIds(); + itk::Point GetItkPoint(double point[3]); + void UpdateFiberGeometry(); + virtual void PrintSelf(std::ostream &os, itk::Indent indent) const override; private: // actual fiber container vtkSmartPointer m_FiberPolyData; // contains fiber ids vtkSmartPointer m_FiberIdDataSet; int m_NumFibers; vtkSmartPointer m_FiberColors; vtkSmartPointer m_FiberWeights; std::vector< float > m_FiberLengths; float m_MinFiberLength; float m_MaxFiberLength; float m_MeanFiberLength; float m_MedianFiberLength; float m_LengthStDev; itk::TimeStamp m_UpdateTime2D; itk::TimeStamp m_UpdateTime3D; mitk::BaseGeometry::Pointer m_ReferenceGeometry; }; } // namespace mitk #endif /* _MITK_FiberBundle_H */ diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp index ecf520ae8f..ef672c393f 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibility.cpp @@ -1,208 +1,217 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; typedef itksys::SystemTools ist; typedef itk::Image ItkUcharImgType; typedef std::tuple< ItkUcharImgType::Pointer, std::string > MaskType; void CreateFolderStructure(const std::string& path) { if (ist::PathExists(path)) ist::RemoveADirectory(path); ist::MakeDirectory(path); ist::MakeDirectory(path + "/known_tracts/"); - ist::MakeDirectory(path + "/plausible_tracts/"); + ist::MakeDirectory(path + "/candidate_tracts/"); ist::MakeDirectory(path + "/implausible_tracts/"); + ist::MakeDirectory(path + "/skipped_masks/"); } ItkUcharImgType::Pointer LoadItkMaskImage(const std::string& filename) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); ItkUcharImgType::Pointer itkMask = ItkUcharImgType::New(); mitk::CastToItkImage(img, itkMask); return itkMask; } -std::vector< MaskType > get_file_list(const std::string& path, const std::string& skipped) +std::vector< MaskType > get_file_list(const std::string& path, int dropout, const std::string& skipped_path) { std::srand(std::time(0)); std::vector< MaskType > mask_list; itk::Directory::Pointer dir = itk::Directory::New(); - ofstream myfile; - myfile.open (skipped); if (dir->Load(path.c_str())) { int n = dir->GetNumberOfFiles(); for (int r = 0; r < n; r++) { const char *filename = dir->GetFile(r); - if (std::rand()%3==0 && ist::GetFilenameWithoutExtension(filename)!="CC") + if (dropout>1 && std::rand()%dropout==0 && ist::GetFilenameWithoutExtension(filename)!="CC") { - MITK_INFO << "Dismissing " << ist::GetFilenameWithoutExtension(filename); - myfile << ist::GetFilenameName(filename) << "\n"; + MITK_INFO << "Skipping " << ist::GetFilenameWithoutExtension(filename); + ist::CopyAFile(path + '/' + filename, skipped_path + ist::GetFilenameName(filename)); continue; } std::string ext = ist::GetFilenameExtension(filename); if (ext==".nii" || ext==".nii.gz" || ext==".nrrd") { MITK_INFO << "Loading " << ist::GetFilenameWithoutExtension(filename); streambuf *old = cout.rdbuf(); // <-- save stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect MaskType m(LoadItkMaskImage(path + '/' + filename), ist::GetFilenameWithoutExtension(filename)); mask_list.push_back(m); std::cout.rdbuf (old); // <-- restore } } } - myfile.close(); return mask_list; } mitk::FiberBundle::Pointer TransformToMRtrixSpace(mitk::FiberBundle::Pointer fib) { 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(fib->GetFiberPolyData()); transformFilter->SetTransform(geometry->GetVtkTransform()); transformFilter->Update(); mitk::FiberBundle::Pointer transformed_fib = mitk::FiberBundle::New(transformFilter->GetOutput()); return transformed_fib; } /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Tract Plausibility"); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::InputFile, "Input:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output folder", us::Any(), false); parser.addArgument("reference_mask_folder", "m", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks of known bundles", false); - parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "", "", false); + parser.addArgument("gray_matter_mask", "gm", mitkCommandLineParser::String, "", ""); + parser.addArgument("dropout", "", mitkCommandLineParser::Int, "", "", false); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string fibFile = us::any_cast(parsedArgs["input"]); - string gray_matter_mask = us::any_cast(parsedArgs["gray_matter_mask"]); string reference_mask_folder = us::any_cast(parsedArgs["reference_mask_folder"]); string out_folder = us::any_cast(parsedArgs["out"]); + string gray_matter_mask = ""; + if (parsedArgs.count("gray_matter_mask")) + gray_matter_mask = us::any_cast(parsedArgs["gray_matter_mask"]); + + int dropout = 1; + if (parsedArgs.count("dropout")) + dropout = us::any_cast(parsedArgs["dropout"]); + try { CreateFolderStructure(out_folder); - std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder, out_folder + "skipped_masks.txt"); + std::vector< MaskType > known_tract_masks = get_file_list(reference_mask_folder, dropout, out_folder + "/skipped_masks/"); if (known_tract_masks.empty()) return EXIT_FAILURE; - ItkUcharImgType::Pointer gm_image = LoadItkMaskImage(gray_matter_mask); mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(fibFile)[0].GetPointer()); - // filter gray matter fibers - mitk::FiberBundle::Pointer not_gm_fibers = inputTractogram->ExtractFiberSubset(gm_image, false, true, true); - mitk::IOUtil::Save(not_gm_fibers, out_folder + "/implausible_tracts/no_gm_endings.trk"); - inputTractogram = inputTractogram->ExtractFiberSubset(gm_image, false, false, true); + if (gray_matter_mask.compare("")!=0) + { + ItkUcharImgType::Pointer gm_image = LoadItkMaskImage(gray_matter_mask); + // filter gray matter fibers + mitk::FiberBundle::Pointer not_gm_fibers = inputTractogram->ExtractFiberSubset(gm_image, false, true, true); + mitk::IOUtil::Save(not_gm_fibers, out_folder + "/implausible_tracts/no_gm_endings.trk"); + inputTractogram = inputTractogram->ExtractFiberSubset(gm_image, false, false, true); + } // resample fibers float minSpacing = 1; if(std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[1] && std::get<0>(known_tract_masks.at(0))->GetSpacing()[0](known_tract_masks.at(0))->GetSpacing()[2]) minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[0]; else if (std::get<0>(known_tract_masks.at(0))->GetSpacing()[1] < std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]) minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[1]; else minSpacing = std::get<0>(known_tract_masks.at(0))->GetSpacing()[2]; inputTractogram->ResampleLinear(minSpacing/5); // find known tracts via overlap match mitk::FiberBundle::Pointer all_known_tracts = nullptr; for ( MaskType mask : known_tract_masks ) { ItkUcharImgType::Pointer mask_image = std::get<0>(mask); std::string mask_name = std::get<1>(mask); mitk::FiberBundle::Pointer known_tract = inputTractogram->ExtractFiberSubset(mask_image, true, false, false, 0.7, false); mitk::IOUtil::Save(known_tract, out_folder + "/known_tracts/" + mask_name + ".trk"); if (all_known_tracts.IsNull()) all_known_tracts = mitk::FiberBundle::New(known_tract->GetFiberPolyData()); else all_known_tracts = all_known_tracts->AddBundle(known_tract); } mitk::IOUtil::Save(all_known_tracts, out_folder + "/known_tracts/all_known_tracts.trk"); mitk::IOUtil::Save(TransformToMRtrixSpace(all_known_tracts), out_folder + "/known_tracts/all_known_tracts_mrtrixspace.fib"); mitk::FiberBundle::Pointer remaining_tracts = inputTractogram->SubtractBundle(all_known_tracts); - mitk::IOUtil::Save(remaining_tracts, out_folder + "/plausible_tracts/remaining_tracts.trk"); - mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), out_folder + "/plausible_tracts/remaining_tracts_mrtrixspace.fib"); + mitk::IOUtil::Save(remaining_tracts, out_folder + "/candidate_tracts/remaining_tracts.trk"); + mitk::IOUtil::Save(TransformToMRtrixSpace(remaining_tracts), out_folder + "/candidate_tracts/remaining_tracts_mrtrixspace.fib"); MITK_INFO << "step_size: " << minSpacing/5; } catch (itk::ExceptionObject e) { std::cout << e; 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/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp index 6b1f1240c0..98b556829c 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/TractographyEvaluation/TractPlausibilityFit.cpp @@ -1,292 +1,453 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace std; typedef itksys::SystemTools ist; typedef itk::Point PointType4; typedef itk::Image< float, 4 > PeakImgType; -std::vector< string > get_file_list(const std::string& path) +std::vector< string > get_file_list(const std::string& path, std::vector< string > extensions={".fib", ".trk"}) { std::vector< string > file_list; itk::Directory::Pointer dir = itk::Directory::New(); if (dir->Load(path.c_str())) { int n = dir->GetNumberOfFiles(); for (int r = 0; r < n; r++) { const char *filename = dir->GetFile(r); std::string ext = ist::GetFilenameExtension(filename); - if (ext==".fib" || ext==".trk") - file_list.push_back(path + '/' + filename); + for (auto e : extensions) + { + if (ext==e) + { + file_list.push_back(path + '/' + filename); + break; + } + } } } return file_list; } +std::vector< mitk::FiberBundle::Pointer > CombineTractograms(std::vector< mitk::FiberBundle::Pointer > reference, std::vector< mitk::FiberBundle::Pointer > candidates, int skip=-1) +{ + std::vector< mitk::FiberBundle::Pointer > fib; + for (auto f : reference) + fib.push_back(f); + + int c = 0; + for (auto f : candidates) + { + if (c!=skip) + fib.push_back(f); + ++c; + } + + return fib; +} + /*! \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[]) { mitkCommandLineParser parser; parser.setTitle(""); parser.setCategory("Fiber Tracking Evaluation"); parser.setDescription(""); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i1", mitkCommandLineParser::InputFile, "Input tractogram:", "input tractogram (.fib, vtk ascii file format)", us::Any(), false); parser.addArgument("", "i2", mitkCommandLineParser::InputFile, "Input peaks:", "input peak image", us::Any(), false); parser.addArgument("", "i3", mitkCommandLineParser::InputFile, "", "", us::Any(), false); - parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any()); + parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output folder:", "output folder", us::Any()); parser.addArgument("mask", "", mitkCommandLineParser::InputFile, "", "", us::Any(), false); parser.addArgument("min_gain", "", mitkCommandLineParser::Float, "Min. gain:", "process stops if remaining bundles don't contribute enough", 0.05); + parser.addArgument("reference_mask_folder", "", mitkCommandLineParser::String, "Reference Mask Folder:", "reference masks for accuracy evaluation", false); + parser.addArgument("subtract", "", mitkCommandLineParser::Bool, "", "", false); + parser.addArgument("normalize_gain", "", mitkCommandLineParser::Bool, "", "", false); parser.addArgument("max_iter", "", mitkCommandLineParser::Int, "Max. iterations:", "maximum number of optimizer iterations", 20); parser.addArgument("bundle_based", "", mitkCommandLineParser::Bool, "Bundle based fit:", "fit one weight per input tractogram/bundle, not for each fiber", false); parser.addArgument("min_g", "", mitkCommandLineParser::Float, "Min. gradient:", "lower termination threshold for gradient magnitude", 1e-5); parser.addArgument("lambda", "", mitkCommandLineParser::Float, "Lambda:", "modifier for regularization", 0.1); parser.addArgument("dont_filter_outliers", "", mitkCommandLineParser::Bool, "Don't filter outliers:", "don't perform second optimization run with an upper weight bound based on the first weight estimation (95% quantile)", false); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; string fib_file = us::any_cast(parsedArgs["i1"]); string peak_file_name = us::any_cast(parsedArgs["i2"]); string candidate_folder = us::any_cast(parsedArgs["i3"]); - string outRoot = us::any_cast(parsedArgs["o"]); + string out_folder = us::any_cast(parsedArgs["o"]); bool single_fib = true; if (parsedArgs.count("bundle_based")) single_fib = !us::any_cast(parsedArgs["bundle_based"]); + bool greedy_add = true; + if (parsedArgs.count("subtract")) + greedy_add = !us::any_cast(parsedArgs["subtract"]); + 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 min_gain = 0.05; if (parsedArgs.count("min_gain")) min_gain = us::any_cast(parsedArgs["min_gain"]); float lambda = 0.1; if (parsedArgs.count("lambda")) lambda = us::any_cast(parsedArgs["lambda"]); + bool normalize_gain = false; + if (parsedArgs.count("normalize_gain")) + normalize_gain = us::any_cast(parsedArgs["normalize_gain"]); + bool filter_outliers = true; if (parsedArgs.count("dont_filter_outliers")) filter_outliers = !us::any_cast(parsedArgs["dont_filter_outliers"]); + string mask_file = ""; if (parsedArgs.count("mask")) mask_file = us::any_cast(parsedArgs["mask"]); + + string reference_mask_folder = ""; + if (parsedArgs.count("reference_mask_folder")) + reference_mask_folder = us::any_cast(parsedArgs["reference_mask_folder"]); + try { itk::TimeProbe clock; clock.Start(); + if (ist::PathExists(out_folder)) + ist::RemoveADirectory(out_folder); + ist::MakeDirectory(out_folder); + + ofstream logfile; + logfile.open (out_folder + "log.txt"); + + itk::ImageFileWriter< PeakImgType >::Pointer peak_image_writer = itk::ImageFileWriter< PeakImgType >::New(); mitk::PreferenceListReaderOptionsFunctor functor = mitk::PreferenceListReaderOptionsFunctor({"Peak Image", "Fiberbundles"}, {}); mitk::Image::Pointer inputImage = dynamic_cast(mitk::IOUtil::Load(peak_file_name, &functor)[0].GetPointer()); + float minSpacing = 1; + if(inputImage->GetGeometry()->GetSpacing()[0]GetGeometry()->GetSpacing()[1] && inputImage->GetGeometry()->GetSpacing()[0]GetGeometry()->GetSpacing()[2]) + minSpacing = inputImage->GetGeometry()->GetSpacing()[0]; + else if (inputImage->GetGeometry()->GetSpacing()[1] < inputImage->GetGeometry()->GetSpacing()[2]) + minSpacing = inputImage->GetGeometry()->GetSpacing()[1]; + else + minSpacing = inputImage->GetGeometry()->GetSpacing()[2]; + + // Load mask file. Fit is only performed inside the mask itk::FitFibersToImageFilter::UcharImgType::Pointer mask = nullptr; if (mask_file.compare("")!=0) { mitk::Image::Pointer mitk_mask = dynamic_cast(mitk::IOUtil::Load(mask_file)[0].GetPointer()); mitk::CastToItkImage(mitk_mask, mask); } + // Load masks covering the true positives for evaluation purposes + std::vector< string > reference_mask_filenames; + std::vector< itk::FitFibersToImageFilter::UcharImgType::Pointer > reference_masks; + if (reference_mask_folder.compare("")!=0) + { + reference_mask_filenames = get_file_list(reference_mask_folder, {".nii", ".nii.gz", ".nrrd"}); + for (auto filename : reference_mask_filenames) + { + itk::FitFibersToImageFilter::UcharImgType::Pointer ref_mask = nullptr; + mitk::Image::Pointer ref_mitk_mask = dynamic_cast(mitk::IOUtil::Load(filename)[0].GetPointer()); + mitk::CastToItkImage(ref_mitk_mask, ref_mask); + reference_masks.push_back(ref_mask); + } + } + + // Load peak image typedef mitk::ImageToItk< PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(inputImage); caster->Update(); PeakImgType::Pointer peak_image = caster->GetOutput(); + // Load reference tractogram consisting of all known tracts std::vector< mitk::FiberBundle::Pointer > input_reference; mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(fib_file)[0].GetPointer()); if (fib.IsNull()) return EXIT_FAILURE; + fib->ResampleLinear(minSpacing/10.0); input_reference.push_back(fib); + // Load all candidate tracts std::vector< mitk::FiberBundle::Pointer > input_candidates; std::vector< string > candidate_tract_files = get_file_list(candidate_folder); for (string f : candidate_tract_files) { mitk::FiberBundle::Pointer fib = dynamic_cast(mitk::IOUtil::Load(f)[0].GetPointer()); if (fib.IsNull()) continue; + fib->ResampleLinear(minSpacing/10.0); input_candidates.push_back(fib); } + // Fit known tracts to peak image to obtain underexplained image int iteration = 0; - std::string name = ist::GetFilenameWithoutExtension(fib_file); - itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); - fitter->SetTractograms(input_reference); + if (greedy_add) + fitter->SetTractograms(input_reference); + else + fitter->SetTractograms(CombineTractograms(input_reference, input_candidates)); fitter->SetFitIndividualFibers(single_fib); fitter->SetMaxIterations(max_iter); fitter->SetGradientTolerance(g_tol); fitter->SetLambda(lambda); fitter->SetFilterOutliers(filter_outliers); fitter->SetPeakImage(peak_image); fitter->SetVerbose(true); - fitter->SetDeepCopy(false); + fitter->SetResampleFibers(false); fitter->SetMaskImage(mask); fitter->Update(); - -// fitter->GetTractograms().at(0)->SetFiberWeights(fitter->GetCoverage()); -// fitter->GetTractograms().at(0)->ColorFibersByFiberWeights(false, false); - - mitk::IOUtil::Save(fitter->GetTractograms().at(0), outRoot + "0_" + name + ".fib"); - peak_image = fitter->GetUnderexplainedImage(); - - itk::ImageFileWriter< PeakImgType >::Pointer writer = itk::ImageFileWriter< PeakImgType >::New(); - writer->SetInput(peak_image); - writer->SetFileName(outRoot + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); - writer->Update(); + std::string name = ist::GetFilenameWithoutExtension(fib_file); + mitk::IOUtil::Save(fitter->GetTractograms().at(0), out_folder + "0_" + name + ".fib"); + if (greedy_add) + { + peak_image = fitter->GetUnderexplainedImage(); + peak_image_writer->SetInput(peak_image); + peak_image_writer->SetFileName(out_folder + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); + peak_image_writer->Update(); + } double coverage = fitter->GetCoverage(); - MITK_INFO << "Iteration: " << iteration; - MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%"; -// fitter->SetPeakImage(peak_image); + MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%"; + logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "%\n\n"; + // fitter->SetPeakImage(peak_image); + // Iteratively add/subtract candidate bundles in a greedy manner while (!input_candidates.empty()) { - streambuf *old = cout.rdbuf(); // <-- save - stringstream ss; - std::cout.rdbuf (ss.rdbuf()); // <-- redirect - double next_coverage = 0; + double next_covered_directions = 1.0; mitk::FiberBundle::Pointer best_candidate = nullptr; - for (auto fib : input_candidates) + if (greedy_add) { - - // WHY NECESSARY AGAIN?? - itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); - fitter->SetFitIndividualFibers(single_fib); - fitter->SetMaxIterations(max_iter); - fitter->SetGradientTolerance(g_tol); - fitter->SetLambda(lambda); - fitter->SetFilterOutliers(filter_outliers); - fitter->SetVerbose(false); - fitter->SetPeakImage(peak_image); - fitter->SetDeepCopy(false); - fitter->SetMaskImage(mask); - // ****************************** - fitter->SetTractograms({fib}); - fitter->Update(); - - double candidate_coverage = fitter->GetCoverage(); - - if (candidate_coverage>next_coverage) + for (auto fib : input_candidates) { - next_coverage = candidate_coverage; - if ((1.0-coverage) * next_coverage >= min_gain) + // WHY NECESSARY AGAIN?? + itk::FitFibersToImageFilter::Pointer fitter = itk::FitFibersToImageFilter::New(); + fitter->SetFitIndividualFibers(single_fib); + fitter->SetMaxIterations(max_iter); + fitter->SetGradientTolerance(g_tol); + fitter->SetLambda(lambda); + fitter->SetFilterOutliers(filter_outliers); + fitter->SetVerbose(false); + fitter->SetPeakImage(peak_image); + fitter->SetResampleFibers(false); + fitter->SetMaskImage(mask); + // ****************************** + fitter->SetTractograms({fib}); + + streambuf *old = cout.rdbuf(); // <-- save + stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + fitter->Update(); + std::cout.rdbuf (old); // <-- restore + + double candidate_coverage = fitter->GetCoverage(); + if (normalize_gain) + candidate_coverage /= fitter->GetNumCoveredDirections(); + + if (candidate_coverage>next_coverage) { - best_candidate = fitter->GetTractograms().at(0); - peak_image = fitter->GetUnderexplainedImage(); + next_covered_directions = fitter->GetNumCoveredDirections(); + next_coverage = candidate_coverage; + if ((1.0-coverage) * next_coverage >= min_gain) + { + best_candidate = fitter->GetTractograms().at(0); + peak_image = fitter->GetUnderexplainedImage(); + } } } } - - if (best_candidate.IsNull()) + else { - std::cout.rdbuf (old); // <-- restore - break; + next_coverage = 1.0; + for (unsigned int i=0; iSetFitIndividualFibers(single_fib); + fitter->SetMaxIterations(max_iter); + fitter->SetGradientTolerance(g_tol); + fitter->SetLambda(lambda); + fitter->SetFilterOutliers(filter_outliers); + fitter->SetVerbose(false); + fitter->SetPeakImage(peak_image); + fitter->SetResampleFibers(false); + fitter->SetMaskImage(mask); + // ****************************** + fitter->SetTractograms(CombineTractograms(input_reference, input_candidates, i)); + + streambuf *old = cout.rdbuf(); // <-- save + stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + fitter->Update(); + std::cout.rdbuf (old); // <-- restore + + double candidate_coverage = fitter->GetCoverage(); + if (normalize_gain) + candidate_coverage /= fitter->GetNumCoveredDirections(); + + if (candidate_coverageGetNumCoveredDirections(); + next_coverage = candidate_coverage; +// if ((1.0-coverage) * next_coverage >= min_gain) +// { + best_candidate = input_candidates.at(i); +// } + } + } } -// fitter->SetPeakImage(peak_image); -// best_candidate->SetFiberWeights((1.0-coverage) * next_coverage); -// best_candidate->ColorFibersByFiberWeights(false, false); + if (best_candidate.IsNull()) + break; - coverage += (1.0-coverage) * next_coverage; + // fitter->SetPeakImage(peak_image); int i=0; std::vector< mitk::FiberBundle::Pointer > remaining_candidates; std::vector< 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++; - mitk::IOUtil::Save(best_candidate, outRoot + boost::lexical_cast(iteration) + "_" + name + ".fib"); - writer->SetInput(peak_image); - writer->SetFileName(outRoot + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); - writer->Update(); + streambuf *old = cout.rdbuf(); // <-- save + stringstream ss; + std::cout.rdbuf (ss.rdbuf()); // <-- redirect + // Save winning candidate + mitk::IOUtil::Save(best_candidate, out_folder + boost::lexical_cast(iteration) + "_" + name + ".fib"); + if (greedy_add) + { + peak_image_writer->SetInput(peak_image); + peak_image_writer->SetFileName(out_folder + boost::lexical_cast(iteration) + "_" + name + ".nrrd"); + peak_image_writer->Update(); + } + + // Calculate best overlap with reference masks for evaluation purposes + float best_overlap = 0; + int best_overlap_index = -1; + i = 0; + for (auto ref_mask : reference_masks) + { + float overlap = best_candidate->GetOverlap(ref_mask, false); + if (overlap>best_overlap) + { + best_overlap = overlap; + best_overlap_index = i; + } + ++i; + } std::cout.rdbuf (old); // <-- restore - MITK_INFO << "Iteration: " << iteration; - MITK_INFO << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "% (+" << 100*(1.0-coverage) * next_coverage << "%)"; + if (normalize_gain) + next_coverage *= next_covered_directions; + if (greedy_add) + { + MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*(coverage + (1.0-coverage) * next_coverage) << "% (+" << 100*(1.0-coverage) * next_coverage << "%)"; + logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*(coverage + (1.0-coverage) * next_coverage) << "% (+" << 100*(1.0-coverage) * next_coverage << "%)\n"; + coverage += (1.0-coverage) * next_coverage; + } + else + { + MITK_INFO << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "% (-" << 100*(coverage-next_coverage) << "%)"; + logfile << "Iteration: " << iteration << " - " << std::fixed << "Coverage: " << setprecision(2) << 100.0*coverage << "% (-" << 100*(coverage-next_coverage) << "%)\n"; + coverage = next_coverage; + } + + if (best_overlap_index>=0) + { + MITK_INFO << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)); + logfile << "Best overlap: " << best_overlap << " - " << ist::GetFilenameWithoutExtension(reference_mask_filenames.at(best_overlap_index)) << "\n"; + } + logfile << "\n"; } clock.Stop(); int h = clock.GetTotal()/3600; int m = ((int)clock.GetTotal()%3600)/60; int s = (int)clock.GetTotal()%60; MITK_INFO << "Plausibility estimation took " << h << "h, " << m << "m and " << s << "s"; + logfile.close(); } catch (itk::ExceptionObject e) { std::cout << e; return EXIT_FAILURE; } catch (std::exception e) { std::cout << e.what(); return EXIT_FAILURE; } catch (...) { std::cout << "ERROR!?!"; return EXIT_FAILURE; } return EXIT_SUCCESS; }