diff --git a/CMake/BuildConfigurations/DiffusionAll.cmake b/CMake/BuildConfigurations/DiffusionAll.cmake index 449fe3c02f..718619ddb0 100644 --- a/CMake/BuildConfigurations/DiffusionAll.cmake +++ b/CMake/BuildConfigurations/DiffusionAll.cmake @@ -1,38 +1,39 @@ 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_MatchPoint ON CACHE BOOL "" FORCE) set(MITK_USE_DCMTK ON CACHE BOOL "" FORCE) set(MITK_USE_DCMQI ON CACHE BOOL "" FORCE) +set(MITK_USE_SWIG ON CACHE BOOL "" FORCE) set(MITK_USE_Python ON CACHE BOOL "" FORCE) set(MITK_USE_BetData ON CACHE BOOL "" FORCE) set(MITK_USE_SYSTEM_PYTHON 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_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/CMake/BuildConfigurations/DiffusionRelease.cmake b/CMake/BuildConfigurations/DiffusionRelease.cmake index 4945db9bc6..b1c58986e7 100644 --- a/CMake/BuildConfigurations/DiffusionRelease.cmake +++ b/CMake/BuildConfigurations/DiffusionRelease.cmake @@ -1,41 +1,42 @@ message(STATUS "Configuring MITK Diffusion Release Build") # 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_MatchPoint ON CACHE BOOL "" FORCE) set(MITK_USE_DCMTK ON CACHE BOOL "" FORCE) set(MITK_USE_DCMQI ON CACHE BOOL "" FORCE) +set(MITK_USE_SWIG ON CACHE BOOL "" FORCE) set(MITK_USE_Python ON CACHE BOOL "" FORCE) set(MITK_USE_BetData ON CACHE BOOL "" FORCE) set(MITK_USE_SYSTEM_PYTHON 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 OFF 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 OFF CACHE BOOL "Build miscellaneous commandline tools for diffusion" FORCE) set(BUILD_DiffusionQuantificationCmdApps OFF 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 OFF 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 OFF 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) # Enable 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) set(CMAKE_BUILD_TYPE Release CACHE STRING "" FORCE) diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp index dfdaf409cb..1abec0e32e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.cpp @@ -1,163 +1,168 @@ /*=================================================================== 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 "itkFiberCurvatureFilter.h" #define _USE_MATH_DEFINES #include #include #include #include namespace itk{ FiberCurvatureFilter::FiberCurvatureFilter() : m_AngularDeviation(30) - , m_Distance(5.0) + , m_Distance(10.0) , m_RemoveFibers(false) + , m_UseMedian(false) { } FiberCurvatureFilter::~FiberCurvatureFilter() { } void FiberCurvatureFilter::GenerateData() { vtkSmartPointer inputPoly = m_InputFiberBundle->GetFiberPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(inputPoly->GetNumberOfCells()); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; #pragma omp critical { ++disp; vtkCell* cell = inputPoly->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p); vnl_vector_fixed< double, 3 > p_vec; p_vec[0]=p[0]; p_vec[1]=p[1]; p_vec[2]=p[2]; vertices.push_back(p_vec); } } // calculate curvatures int numPoints = vertices.size(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; j > vectors; vnl_vector_fixed< double, 3 > meanV; meanV.fill(0.0); while(dist1) { vnl_vector_fixed< double, 3 > p1 = vertices.at(c-1); vnl_vector_fixed< double, 3 > p2 = vertices.at(c); vnl_vector_fixed< double, 3 > v = p2-p1; dist += v.magnitude(); v.normalize(); vectors.push_back(v); - if (c==j) + if (m_UseMedian && c==j) meanV += v; + else if (!m_UseMedian) + meanV += v; c--; } c = j; dist = 0; while(dist p1 = vertices.at(c); vnl_vector_fixed< double, 3 > p2 = vertices.at(c+1); vnl_vector_fixed< double, 3 > v = p2-p1; dist += v.magnitude(); v.normalize(); vectors.push_back(v); - if (c==j) + if (m_UseMedian && c==j) meanV += v; + else if (!m_UseMedian) + meanV += v; c++; } meanV.normalize(); double dev = 0; for (auto vec : vectors) { double angle = dot_product(meanV, vec); if (angle>1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); if (devInsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } else { if (m_RemoveFibers) { container = vtkSmartPointer::New(); break; } if (container->GetNumberOfPoints()>0) { #pragma omp critical vtkNewCells->InsertNextCell(container); } container = vtkSmartPointer::New(); } } #pragma omp critical { if (container->GetNumberOfPoints()>0) vtkNewCells->InsertNextCell(container); } } vtkSmartPointer outputPoly = vtkSmartPointer::New(); outputPoly->SetPoints(vtkNewPoints); outputPoly->SetLines(vtkNewCells); m_OutputFiberBundle = FiberBundle::New(outputPoly); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h index 82c84a06a4..effb674dd4 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkFiberCurvatureFilter.h @@ -1,81 +1,83 @@ /*=================================================================== 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 itkFiberCurvatureFilter_h #define itkFiberCurvatureFilter_h // MITK #include #include #include // ITK #include // VTK #include #include #include #include #include namespace itk{ /** * \brief */ class FiberCurvatureFilter : public ProcessObject { public: typedef FiberCurvatureFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkTypeMacro( FiberCurvatureFilter, ProcessObject ) virtual void Update() override{ this->GenerateData(); } itkSetMacro( Distance, double ) itkSetMacro( AngularDeviation, double ) itkSetMacro( RemoveFibers, bool ) + itkSetMacro( UseMedian, bool ) itkSetMacro( InputFiberBundle, FiberBundle::Pointer ) itkGetMacro( OutputFiberBundle, FiberBundle::Pointer ) protected: void GenerateData() override; FiberCurvatureFilter(); virtual ~FiberCurvatureFilter(); FiberBundle::Pointer m_InputFiberBundle; FiberBundle::Pointer m_OutputFiberBundle; double m_AngularDeviation; double m_Distance; bool m_RemoveFibers; + bool m_UseMedian; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkFiberCurvatureFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 1e5c3dd017..6510d699ef 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,906 +1,890 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include "itkStreamlineTrackingFilter.h" #include #include #include #include #include "itkPointShell.h" #include #include +#include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_BuildFibersFinished(false) , m_BuildFibersReady(0) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_StoppingRegions(nullptr) , m_TargetRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_OutputProbabilityMap(nullptr) + , m_MinVoxelSize(-1) , m_AngularThresholdDeg(-1) , m_StepSizeVox(-1) , m_SamplingDistanceVox(-1) , m_AngularThreshold(-1) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) , m_AvoidStop(true) , m_RandomSampling(false) , m_SamplingDistance(-1) , m_DeflectionMod(1.0) , m_OnlyForwardSamples(true) , m_UseStopVotes(true) , m_NumberOfSamples(30) , m_NumPreviousDirections(1) , m_MaxNumTracts(-1) , m_Verbose(true) - , m_AposterioriCurvCheck(false) + , m_LoopCheck(-1) , m_DemoMode(false) , m_Random(true) , m_UseOutputProbabilityMap(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) , m_InterpolateMask(true) + , m_TrialsPerSeed(10) { this->SetNumberOfRequiredInputs(0); } std::string StreamlineTrackingFilter::GetStatusText() { std::string status = "Seedpoints processed: " + boost::lexical_cast(m_Progress) + "/" + boost::lexical_cast(m_SeedPoints.size()); if (m_SeedPoints.size()>0) status += " (" + boost::lexical_cast(100*m_Progress/m_SeedPoints.size()) + "%)"; if (m_MaxNumTracts>0) status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts) + "/" + boost::lexical_cast(m_MaxNumTracts); else status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts); return status; } void StreamlineTrackingFilter::BeforeTracking() { m_StopTracking = false; m_TrackingHandler->SetRandom(m_Random); m_TrackingHandler->InitForTracking(); m_FiberPolyData = PolyDataType::New(); m_Points = vtkSmartPointer< vtkPoints >::New(); m_Cells = vtkSmartPointer< vtkCellArray >::New(); itk::Vector< double, 3 > imageSpacing = m_TrackingHandler->GetSpacing(); - float minSpacing; if(imageSpacing[0]SetAngularThreshold(m_AngularThreshold); if (m_SamplingDistanceVoxGetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } if (m_UseOutputProbabilityMap) { m_OutputProbabilityMap = ItkDoubleImgType::New(); m_OutputProbabilityMap->SetSpacing(imageSpacing); m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin()); m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection()); m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion()); m_OutputProbabilityMap->Allocate(); m_OutputProbabilityMap->FillBuffer(0); } m_MaskInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_StopInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_SeedInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); m_TargetInterpolator = itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::New(); if (m_StoppingRegions.IsNull()) { m_StoppingRegions = ItkFloatImgType::New(); m_StoppingRegions->SetSpacing( imageSpacing ); m_StoppingRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_StoppingRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_StoppingRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_StoppingRegions->Allocate(); m_StoppingRegions->FillBuffer(0); } else std::cout << "StreamlineTracking - Using stopping region image" << std::endl; m_StopInterpolator->SetInputImage(m_StoppingRegions); if (m_TargetRegions.IsNull()) { m_TargetImageSet = false; m_TargetRegions = ItkFloatImgType::New(); m_TargetRegions->SetSpacing( imageSpacing ); m_TargetRegions->SetOrigin( m_TrackingHandler->GetOrigin() ); m_TargetRegions->SetDirection( m_TrackingHandler->GetDirection() ); m_TargetRegions->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_TargetRegions->Allocate(); m_TargetRegions->FillBuffer(1); } else { m_TargetImageSet = true; m_TargetInterpolator->SetInputImage(m_TargetRegions); std::cout << "StreamlineTracking - Using target region image" << std::endl; } if (m_SeedImage.IsNull()) { m_SeedImageSet = false; m_SeedImage = ItkFloatImgType::New(); m_SeedImage->SetSpacing( imageSpacing ); m_SeedImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_SeedImage->SetDirection( m_TrackingHandler->GetDirection() ); m_SeedImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_SeedImage->Allocate(); m_SeedImage->FillBuffer(1); } else { m_SeedImageSet = true; std::cout << "StreamlineTracking - Using seed image" << std::endl; } m_SeedInterpolator->SetInputImage(m_SeedImage); if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkFloatImgType::New(); m_MaskImage->SetSpacing( imageSpacing ); m_MaskImage->SetOrigin( m_TrackingHandler->GetOrigin() ); m_MaskImage->SetDirection( m_TrackingHandler->GetDirection() ); m_MaskImage->SetRegions( m_TrackingHandler->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } else std::cout << "StreamlineTracking - Using mask image" << std::endl; m_MaskInterpolator->SetInputImage(m_MaskImage); if (m_SeedPoints.empty()) GetSeedPointsFromSeedImage(); m_BuildFibersReady = 0; m_BuildFibersFinished = false; m_Tractogram.clear(); m_SamplingPointset = mitk::PointSet::New(); m_AlternativePointset = mitk::PointSet::New(); m_StopVotePointset = mitk::PointSet::New(); m_StartTime = std::chrono::system_clock::now(); if (m_DemoMode) omp_set_num_threads(1); if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::DETERMINISTIC) std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; else if(m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::PROBABILISTIC) + { std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; + std::cout << "StreamlineTracking - Trials per seed: " << m_TrialsPerSeed << std::endl; + } else std::cout << "StreamlineTracking - Mode: ???" << std::endl; std::cout << "StreamlineTracking - Angular threshold: " << m_AngularThreshold << " (" << 180*std::acos( m_AngularThreshold )/M_PI << "°)" << std::endl; - std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/minSpacing << "*vox)" << std::endl; + std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/m_MinVoxelSize << "*vox)" << std::endl; std::cout << "StreamlineTracking - Seeds per voxel: " << m_SeedsPerVoxel << std::endl; std::cout << "StreamlineTracking - Max. tract length: " << m_MaxTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Min. tract length: " << m_MinTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Max. num. tracts: " << m_MaxNumTracts << std::endl; + std::cout << "StreamlineTracking - Loop check: " << m_LoopCheck << "°" << std::endl; std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_NumberOfSamples << std::endl; - std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/minSpacing << "*vox)" << std::endl; + std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/m_MinVoxelSize << "*vox)" << std::endl; std::cout << "StreamlineTracking - Deflection modifier: " << m_DeflectionMod << std::endl; std::cout << "StreamlineTracking - Use stop votes: " << m_UseStopVotes << std::endl; std::cout << "StreamlineTracking - Only frontal samples: " << m_OnlyForwardSamples << std::endl; if (m_DemoMode) { std::cout << "StreamlineTracking - Running in demo mode"; std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl; } else std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl; } void StreamlineTrackingFilter::CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir) { pos[0] += dir[0]*m_StepSize; pos[1] += dir[1]*m_StepSize; pos[2] += dir[2]*m_StepSize; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< float > theta; theta.resize(NPoints); std::vector< float > phi; phi.resize(NPoints); float C = sqrt(4*M_PI); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(int i=0; i0 && i d; d[0] = cos(theta[i]) * cos(phi[i]); d[1] = cos(theta[i]) * sin(phi[i]); d[2] = sin(theta[i]); pointshell.push_back(d); } return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { if (m_DemoMode) { m_SamplingPointset->Clear(); m_AlternativePointset->Clear(); m_StopVotePointset->Clear(); } vnl_vector_fixed direction; direction.fill(0); if (mitk::imv::IsInsideMask(pos, m_InterpolateMask, m_MaskInterpolator) && !mitk::imv::IsInsideMask(pos, m_InterpolateMask, m_StopInterpolator)) direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position else return direction; vnl_vector_fixed olddir = olddirs.back(); std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_NumberOfSamples); itk::Point sample_pos; int alternatives = 1; int stop_votes = 0; int possible_stop_votes = 0; for (unsigned int i=0; i d; bool is_stop_voter = false; if (m_Random && m_RandomSampling) { d[0] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[1] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[2] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d.normalize(); d *= m_TrackingHandler->GetRandDouble(0,m_SamplingDistance); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); if (m_UseStopVotes && dot>0.7) { is_stop_voter = true; possible_stop_votes++; } else if (m_OnlyForwardSamples && dot<0) continue; d *= m_SamplingDistance; } sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (mitk::imv::IsInsideMask(sample_pos, m_InterpolateMask, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) { direction += tempDir; if(m_DemoMode) m_SamplingPointset->InsertPoint(i, sample_pos); } else if (m_AvoidStop && olddir.magnitude()>0.5) // out of white matter { if (is_stop_voter) stop_votes++; if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); float dot = dot_product(d, olddir); if (dot >= 0.0) // in front of plane defined by pos and olddir d = -d + 2*dot*olddir; // reflect else d = -d; // invert // look a bit further into the other direction sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; alternatives++; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (mitk::imv::IsInsideMask(sample_pos, m_InterpolateMask, m_MaskInterpolator)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) // are we back in the white matter? { direction += d * m_DeflectionMod; // go into the direction of the white matter direction += tempDir; // go into the direction of the white matter direction at this location if(m_DemoMode) m_AlternativePointset->InsertPoint(alternatives, sample_pos); } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); } } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); if (is_stop_voter) stop_votes++; } } if (direction.magnitude()>0.001 && (possible_stop_votes==0 || (float)stop_votes/possible_stop_votes<0.5) ) direction.normalize(); else direction.fill(0); return direction; } -float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front) +float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, DirectionContainer* container, float tractLength, bool front) { vnl_vector_fixed zero_dir; zero_dir.fill(0.0); std::deque< vnl_vector_fixed > last_dirs; for (unsigned int i=0; i oldIndex; m_TrackingHandler->WorldToIndex(pos, oldIndex); // get new position CalculateNewPosition(pos, dir); // is new position inside of image and mask if (m_AbortTracking) // if not end streamline - { return tractLength; + + // if yes, add new point to streamline + dir.normalize(); + if (front) + { + fib->push_front(pos); + container->push_front(dir); } - else // if yes, add new point to streamline + else { - tractLength += m_StepSize; - if (front) - fib->push_front(pos); - else - fib->push_back(pos); + fib->push_back(pos); + container->push_back(dir); + } + tractLength += m_StepSize; - if (m_AposterioriCurvCheck) - { - int curv = CheckCurvature(fib, front); - if (curv>0) - { - tractLength -= m_StepSize*curv; - while (curv>0) - { - if (front) - fib->pop_front(); - else - fib->pop_back(); - curv--; - } - return tractLength; - } - } + if (m_LoopCheck>=0 && CheckCurvature(container, front)>m_LoopCheck) + return tractLength; - if (tractLength>m_MaxTractLength) - return tractLength; - } + if (tractLength>m_MaxTractLength) + return tractLength; if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? { #pragma omp critical { m_BuildFibersReady++; m_Tractogram.push_back(*fib); BuildFibers(true); m_Stop = true; while (m_Stop){ } } } - dir.normalize(); last_dirs.push_back(dir); if (last_dirs.size()>m_NumPreviousDirections) last_dirs.pop_front(); dir = GetNewDirection(pos, last_dirs, oldIndex); while (m_PauseTracking){} if (dir.magnitude()<0.0001) return tractLength; } return tractLength; } -int StreamlineTrackingFilter::CheckCurvature(FiberType* fib, bool front) +float StreamlineTrackingFilter::CheckCurvature(DirectionContainer* fib, bool front) { - float m_Distance = 5; if (fib->size()<3) return 0; - + float m_Distance = m_MinVoxelSize*4; float dist = 0; + std::vector< vnl_vector_fixed< float, 3 > > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); float dev = 0; if (front) { - int c=0; + int c = 0; while(distsize()-1) { - itk::Point p1 = fib->at(c); - itk::Point p2 = fib->at(c+1); - - 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(); + dist += m_StepSize; + vnl_vector_fixed< float, 3 > v = fib->at(c); vectors.push_back(v); - if (c==0) - meanV += v; + meanV += v; c++; } } else { - int c=fib->size()-1; - while(dist0) + int c = fib->size()-1; + while(dist=0) { - itk::Point p1 = fib->at(c); - itk::Point p2 = fib->at(c-1); - - 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(); + dist += m_StepSize; + vnl_vector_fixed< float, 3 > v = fib->at(c); vectors.push_back(v); - if (c==(int)fib->size()-1) - meanV += v; + meanV += v; c--; } } meanV.normalize(); 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(); - if (dev<30) - return 0; - else - return vectors.size(); + return dev; } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - Calculating seed points."; m_SeedPoints.clear(); typedef ImageRegionConstIterator< ItkFloatImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); sit.GoToBegin(); while (!sit.IsAtEnd()) { if (sit.Value()>0) { ItkFloatImgType::IndexType index = sit.GetIndex(); itk::ContinuousIndex start; start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); if ( mitk::imv::IsInsideMask(worldPos, m_InterpolateMask, m_MaskInterpolator) ) { m_SeedPoints.push_back(worldPos); for (int s = 1; s < m_SeedsPerVoxel; s++) { start[0] = index[0] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[1] = index[1] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[2] = index[2] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); m_SeedPoints.push_back(worldPos); } } } ++sit; } } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (m_Random) std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); m_CurrentTracts = 0; int num_seeds = m_SeedPoints.size(); itk::Index<3> zeroIndex; zeroIndex.Fill(0); m_Progress = 0; int i = 0; int print_interval = num_seeds/100; if (print_interval<100) m_Verbose=false; #pragma omp parallel while (i=num_seeds || m_StopTracking) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { m_Progress += print_interval; std::cout << " \r"; if (m_MaxNumTracts>0) std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_MaxNumTracts << '\r'; else std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << '\r'; cout.flush(); } const itk::Point worldPos = m_SeedPoints.at(temp_i); - FiberType fib; - float tractLength = 0; - unsigned int counter = 0; - - // get starting direction - vnl_vector_fixed dir; dir.fill(0.0); - std::deque< vnl_vector_fixed > olddirs; - while (olddirs.size() gm_start_dir; - // if (m_ControlGmEndings) - // { - // gm_start_dir[0] = m_GmStubs[temp_i][1][0] - m_GmStubs[temp_i][0][0]; - // gm_start_dir[1] = m_GmStubs[temp_i][1][1] - m_GmStubs[temp_i][0][1]; - // gm_start_dir[2] = m_GmStubs[temp_i][1][2] - m_GmStubs[temp_i][0][2]; - // gm_start_dir.normalize(); - // olddirs.pop_back(); - // olddirs.push_back(gm_start_dir); - // } - - if (mitk::imv::IsInsideMask(worldPos, m_InterpolateMask, m_MaskInterpolator)) - dir = m_TrackingHandler->ProposeDirection(worldPos, olddirs, zeroIndex); - - if (dir.magnitude()>0.0001) + + for (unsigned int trials=0; trials dir; dir.fill(0.0); + std::deque< vnl_vector_fixed > olddirs; + while (olddirs.size() gm_start_dir; + // if (m_ControlGmEndings) + // { + // gm_start_dir[0] = m_GmStubs[temp_i][1][0] - m_GmStubs[temp_i][0][0]; + // gm_start_dir[1] = m_GmStubs[temp_i][1][1] - m_GmStubs[temp_i][0][1]; + // gm_start_dir[2] = m_GmStubs[temp_i][1][2] - m_GmStubs[temp_i][0][2]; + // gm_start_dir.normalize(); + // olddirs.pop_back(); + // olddirs.push_back(gm_start_dir); + // } + + if (mitk::imv::IsInsideMask(worldPos, m_InterpolateMask, m_MaskInterpolator)) + dir = m_TrackingHandler->ProposeDirection(worldPos, olddirs, zeroIndex); + + bool success = false; + if (dir.magnitude()>0.0001) + { + /// START DIR + // if (m_ControlGmEndings) + // { + // float a = dot_product(gm_start_dir, dir); + // if (a<0) + // dir = -dir; + // } - // forward tracking - tractLength = FollowStreamline(worldPos, dir, &fib, 0, false); - fib.push_front(worldPos); + // forward tracking + tractLength = FollowStreamline(worldPos, dir, &fib, &direction_container, 0, false); + fib.push_front(worldPos); - // backward tracking (only if we don't explicitely start in the GM) - tractLength = FollowStreamline(worldPos, -dir, &fib, tractLength, true); + // backward tracking (only if we don't explicitely start in the GM) + tractLength = FollowStreamline(worldPos, -dir, &fib, &direction_container, tractLength, true); - counter = fib.size(); + counter = fib.size(); - if (tractLength>=m_MinTractLength && counter>=2) - { -#pragma omp critical - if ( IsValidFiber(&fib) ) + if (tractLength>=m_MinTractLength && counter>=2) { - if (!m_StopTracking) - { - if (!m_UseOutputProbabilityMap) - m_Tractogram.push_back(fib); - else - FiberToProbmap(&fib); - m_CurrentTracts++; - } - if (m_MaxNumTracts > 0 && m_CurrentTracts>=static_cast(m_MaxNumTracts)) +#pragma omp critical + if ( IsValidFiber(&fib) ) { if (!m_StopTracking) { - std::cout << " \r"; - MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; + if (!m_UseOutputProbabilityMap) + m_Tractogram.push_back(fib); + else + FiberToProbmap(&fib); + m_CurrentTracts++; + success = true; + } + if (m_MaxNumTracts > 0 && m_CurrentTracts>=static_cast(m_MaxNumTracts)) + { + if (!m_StopTracking) + { + std::cout << " \r"; + MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; + } + m_StopTracking = true; } - m_StopTracking = true; } } } - } - } + if (success || m_TrackingHandler->GetMode()!=mitk::TrackingDataHandler::PROBABILISTIC) + break; // we only try one seed point multiple times if we use a probabilistic tracker and have not found a valid streamline yet + + }// trials per seed + + }// seed points this->AfterTracking(); } bool StreamlineTrackingFilter::IsValidFiber(FiberType* fib) { if (m_TargetImageSet && m_SeedImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMask, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMask, m_TargetInterpolator) ) return true; if ( mitk::imv::IsInsideMask(fib->back(), m_InterpolateMask, m_SeedInterpolator) && mitk::imv::IsInsideMask(fib->front(), m_InterpolateMask, m_TargetInterpolator) ) return true; return false; } else if (m_TargetImageSet) { if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMask, m_TargetInterpolator) && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMask, m_TargetInterpolator) ) return true; return false; } return true; } void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib) { ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0); for (auto p : *fib) { ItkDoubleImgType::IndexType idx; m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); if (idx != last_idx) { if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(idx)) m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1); last_idx = idx; } } } void StreamlineTrackingFilter::BuildFibers(bool check) { if (m_BuildFibersReady::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); for (unsigned int i=0; i container = vtkSmartPointer::New(); FiberType fib = m_Tractogram.at(i); for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) { vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if (check) for (int i=0; iSetPoints(vNewPoints); m_FiberPolyData->SetLines(vNewLines); m_BuildFibersFinished = true; } void StreamlineTrackingFilter::AfterTracking() { if (m_Verbose) std::cout << " \r"; if (!m_UseOutputProbabilityMap) { MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; MITK_INFO << "Generating polydata "; BuildFibers(false); } else { itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer filter = itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); filter->SetInput(m_OutputProbabilityMap); filter->SetOutputMaximum(1.0); filter->SetOutputMinimum(0.0); filter->Update(); m_OutputProbabilityMap = filter->GetOutput(); } MITK_INFO << "done"; m_EndTime = std::chrono::system_clock::now(); std::chrono::hours hh = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::minutes mm = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::seconds ss = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; ss %= 60; MITK_INFO << "Tracking took " << hh.count() << "h, " << mm.count() << "m and " << ss.count() << "s"; m_SeedPoints.clear(); } void StreamlineTrackingFilter::SetDicomProperties(mitk::FiberBundle::Pointer fib) { std::string model_code_value = "-"; std::string model_code_meaning = "-"; std::string algo_code_value = "-"; std::string algo_code_meaning = "-"; if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC && dynamic_cast(m_TrackingHandler) && !m_TrackingHandler->GetInterpolate()) { algo_code_value = "sup181_ee04"; algo_code_meaning = "FACT"; } else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::DETERMINISTIC) { algo_code_value = "sup181_ee01"; algo_code_meaning = "Deterministic"; } else if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::PROBABILISTIC) { algo_code_value = "sup181_ee02"; algo_code_meaning = "Probabilistic"; } if (dynamic_cast(m_TrackingHandler) || (dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetIsOdfFromTensor() ) ) { if ( dynamic_cast(m_TrackingHandler) && dynamic_cast(m_TrackingHandler)->GetNumTensorImages()>1 ) { model_code_value = "sup181_bb02"; model_code_meaning = "Multi Tensor"; } else { model_code_value = "sup181_bb01"; model_code_meaning = "Single Tensor"; } } else if (dynamic_cast*>(m_TrackingHandler) || dynamic_cast*>(m_TrackingHandler)) { model_code_value = "sup181_bb03"; model_code_meaning = "Model Free"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "ODF"; } else if (dynamic_cast(m_TrackingHandler)) { model_code_value = "-"; model_code_meaning = "Peaks"; } fib->SetProperty("DICOM.anatomy.value", mitk::StringProperty::New("T-A0095")); fib->SetProperty("DICOM.anatomy.meaning", mitk::StringProperty::New("White matter of brain and spinal cord")); fib->SetProperty("DICOM.algo_code.value", mitk::StringProperty::New(algo_code_value)); fib->SetProperty("DICOM.algo_code.meaning", mitk::StringProperty::New(algo_code_meaning)); fib->SetProperty("DICOM.model_code.value", mitk::StringProperty::New(model_code_value)); fib->SetProperty("DICOM.model_code.meaning", mitk::StringProperty::New(model_code_meaning)); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h index 8fc5b836a2..d545a62137 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,224 +1,228 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkMLBSTrackingFilter_h_ #define __itkMLBSTrackingFilter_h_ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include namespace itk{ /** * \brief Performs streamline tracking on the input image. Depending on the tracking handler this can be a tensor, peak or machine learning based tracking. */ class MITKFIBERTRACKING_EXPORT StreamlineTrackingFilter : public ProcessObject { public: typedef StreamlineTrackingFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ProcessObject Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(MLBSTrackingFilter, ImageToImageFilter) typedef itk::Image ItkUcharImgType; typedef itk::Image ItkUintImgType; typedef itk::Image ItkDoubleImgType; typedef itk::Image ItkFloatImgType; typedef vtkSmartPointer< vtkPolyData > PolyDataType; + typedef std::deque< vnl_vector_fixed > DirectionContainer; typedef std::deque< itk::Point > FiberType; typedef std::vector< FiberType > BundleType; volatile bool m_PauseTracking; bool m_AbortTracking; bool m_BuildFibersFinished; int m_BuildFibersReady; volatile bool m_Stop; mitk::PointSet::Pointer m_SamplingPointset; mitk::PointSet::Pointer m_StopVotePointset; mitk::PointSet::Pointer m_AlternativePointset; void SetStepSize(float v) ///< Integration step size in voxels, default is 0.5 * voxel { m_StepSizeVox = v; } void SetAngularThreshold(float v) ///< Angular threshold per step (in degree), default is 90deg x stepsize { m_AngularThresholdDeg = v; } void SetSamplingDistance(float v) ///< Maximum distance of sampling points in voxels, default is 0.25 * voxel { m_SamplingDistanceVox = v; } void SetDicomProperties(mitk::FiberBundle::Pointer fib); itkGetMacro( OutputProbabilityMap, ItkDoubleImgType::Pointer) ///< Output probability map itkGetMacro( FiberPolyData, PolyDataType ) ///< Output fibers itkGetMacro( UseOutputProbabilityMap, bool) itkSetMacro( SeedImage, ItkFloatImgType::Pointer) ///< Seeds are only placed inside of this mask. itkSetMacro( MaskImage, ItkFloatImgType::Pointer) ///< Tracking is only performed inside of this mask image. itkSetMacro( SeedsPerVoxel, int) ///< One seed placed in the center of each voxel or multiple seeds randomly placed inside each voxel. itkSetMacro( MinTractLength, float ) ///< Shorter tracts are discarded. itkSetMacro( MaxTractLength, float ) ///< Streamline progression stops if tract is longer than specified. itkSetMacro( UseStopVotes, bool ) ///< Frontal sampling points can vote for stopping the streamline even if the remaining sampling points keep pushing itkSetMacro( OnlyForwardSamples, bool ) ///< Don't use sampling points behind the current position in progression direction itkSetMacro( DeflectionMod, float ) ///< Deflection distance modifier itkSetMacro( StoppingRegions, ItkFloatImgType::Pointer) ///< Streamlines entering a stopping region will stop immediately itkSetMacro( TargetRegions, ItkFloatImgType::Pointer) ///< Only streamline starting and ending in this mask are retained itkSetMacro( DemoMode, bool ) itkSetMacro( NumberOfSamples, unsigned int ) ///< Number of neighborhood sampling points - itkSetMacro( AposterioriCurvCheck, bool ) ///< Checks fiber curvature (angular deviation across 5mm) is larger than 30°. If yes, the streamline progression is stopped. + itkSetMacro( LoopCheck, float ) ///< Checks fiber curvature (angular deviation across 5mm) is larger than 30°. If yes, the streamline progression is stopped. itkSetMacro( AvoidStop, bool ) ///< Use additional sampling points to avoid premature streamline termination itkSetMacro( RandomSampling, bool ) ///< If true, the sampling points are distributed randomly around the current position, not sphericall in the specified sampling distance. itkSetMacro( NumPreviousDirections, unsigned int ) ///< How many "old" steps do we want to consider in our decision where to go next? itkSetMacro( MaxNumTracts, int ) ///< Tracking is stopped if the maximum number of tracts is exceeded itkSetMacro( Random, bool ) ///< If true, seedpoints are shuffled randomly before tracking itkSetMacro( Verbose, bool ) ///< If true, output tracking progress (might be slower) itkSetMacro( UseOutputProbabilityMap, bool) ///< If true, no tractogram but a probability map is created as output. itkSetMacro( StopTracking, bool ) itkSetMacro( InterpolateMask, bool ) + itkSetMacro( TrialsPerSeed, unsigned int ) ///< When using probabilistic tractography, each seed point is used N times until a valid streamline that is compliant with all thresholds etc. is found ///< Use manually defined points in physical space as seed points instead of seed image void SetSeedPoints( const std::vector< itk::Point >& sP) { m_SeedPoints = sP; } void SetTrackingHandler( mitk::TrackingDataHandler* h ) ///< { m_TrackingHandler = h; } virtual void Update() override{ this->GenerateData(); } std::string GetStatusText(); protected: void GenerateData() override; StreamlineTrackingFilter(); ~StreamlineTrackingFilter() {} bool IsValidFiber(FiberType* fib); ///< Check endpoints void FiberToProbmap(FiberType* fib); void GetSeedPointsFromSeedImage(); void CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir); ///< Calculate next integration step. - float FollowStreamline(itk::Point start_pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front); ///< Start streamline in one direction. + float FollowStreamline(itk::Point start_pos, vnl_vector_fixed dir, FiberType* fib, DirectionContainer* container, float tractLength, bool front); ///< Start streamline in one direction. vnl_vector_fixed GetNewDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< Determine new direction by sample voting at the current position taking the last progression direction into account. std::vector< vnl_vector_fixed > CreateDirections(int NPoints); void BeforeTracking(); void AfterTracking(); PolyDataType m_FiberPolyData; vtkSmartPointer m_Points; vtkSmartPointer m_Cells; BundleType m_Tractogram; BundleType m_GmStubs; ItkFloatImgType::Pointer m_StoppingRegions; ItkFloatImgType::Pointer m_TargetRegions; ItkFloatImgType::Pointer m_SeedImage; ItkFloatImgType::Pointer m_MaskImage; ItkDoubleImgType::Pointer m_OutputProbabilityMap; + float m_MinVoxelSize; float m_AngularThresholdDeg; float m_StepSizeVox; float m_SamplingDistanceVox; float m_AngularThreshold; float m_StepSize; int m_MaxLength; float m_MinTractLength; float m_MaxTractLength; int m_SeedsPerVoxel; bool m_AvoidStop; bool m_RandomSampling; float m_SamplingDistance; float m_DeflectionMod; bool m_OnlyForwardSamples; bool m_UseStopVotes; unsigned int m_NumberOfSamples; unsigned int m_NumPreviousDirections; int m_MaxNumTracts; bool m_Verbose; - bool m_AposterioriCurvCheck; + float m_LoopCheck; bool m_DemoMode; bool m_Random; bool m_UseOutputProbabilityMap; std::vector< itk::Point > m_SeedPoints; unsigned int m_CurrentTracts; unsigned int m_Progress; bool m_StopTracking; bool m_InterpolateMask; + unsigned int m_TrialsPerSeed; void BuildFibers(bool check); - int CheckCurvature(FiberType* fib, bool front); + float CheckCurvature(DirectionContainer *fib, bool front); // decision forest mitk::TrackingDataHandler* m_TrackingHandler; std::vector< PolyDataType > m_PolyDataContainer; std::chrono::time_point m_StartTime; std::chrono::time_point m_EndTime; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_MaskInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_StopInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_TargetInterpolator; itk::LinearInterpolateImageFunction< ItkFloatImgType, float >::Pointer m_SeedInterpolator; bool m_SeedImageSet; bool m_TargetImageSet; private: }; } //#ifndef ITK_MANUAL_INSTANTIATION //#include "itkMLBSTrackingFilter.cpp" //#endif #endif //__itkMLBSTrackingFilter_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp index c7feec9017..38cb17d5da 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkTractClusteringFilter.cpp @@ -1,492 +1,497 @@ /*=================================================================== 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 "itkTractClusteringFilter.h" #define _USE_MATH_DEFINES #include #include #include namespace itk{ TractClusteringFilter::TractClusteringFilter() : m_NumPoints(12) , m_InCentroids(nullptr) , m_MinClusterSize(1) , m_MaxClusters(0) , m_MergeDuplicateThreshold(-1) , m_DoResampling(true) , m_FilterMask(nullptr) , m_OverlapThreshold(0.0) { } TractClusteringFilter::~TractClusteringFilter() { for (auto m : m_Metrics) delete m; } std::vector > TractClusteringFilter::GetOutFiberIndices() const { return m_OutFiberIndices; } void TractClusteringFilter::SetMetrics(const std::vector &Metrics) { m_Metrics = Metrics; } std::vector TractClusteringFilter::GetOutClusters() const { return m_OutClusters; } std::vector TractClusteringFilter::GetOutCentroids() const { return m_OutCentroids; } std::vector TractClusteringFilter::GetOutTractograms() const { return m_OutTractograms; } void TractClusteringFilter::SetDistances(const std::vector &Distances) { m_Distances = Distances; } float TractClusteringFilter::CalcOverlap(vnl_matrix& t) { float overlap = 0; if (m_FilterMask.IsNotNull()) { for (unsigned int i=0; i p = t.get_column(i); itk::Point point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; itk::Index<3> idx; m_FilterMask->TransformPhysicalPointToIndex(point, idx); if (m_FilterMask->GetLargestPossibleRegion().IsInside(idx) && m_FilterMask->GetPixel(idx)>0) overlap += 1; } overlap /= m_NumPoints; } else return 1.0; return overlap; } std::vector > TractClusteringFilter::ResampleFibers(mitk::FiberBundle::Pointer tractogram) { mitk::FiberBundle::Pointer temp_fib = tractogram->GetDeepCopy(); if (m_DoResampling) temp_fib->ResampleToNumPoints(m_NumPoints); std::vector< vnl_matrix > out_fib; for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = temp_fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vnl_matrix streamline; streamline.set_size(3, m_NumPoints); streamline.fill(0.0); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< float, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; streamline.set_column(j, candV); } out_fib.push_back(streamline); } return out_fib; } std::vector< TractClusteringFilter::Cluster > TractClusteringFilter::ClusterStep(std::vector< long > f_indices, std::vector distances) { float dist_thres = distances.back(); distances.pop_back(); std::vector< Cluster > C; int N = f_indices.size(); Cluster c1; c1.I.push_back(f_indices.at(0)); c1.h = T[f_indices.at(0)]; c1.n = 1; C.push_back(c1); if (f_indices.size()==1) return C; for (int i=1; i t = T.at(f_indices.at(i)); int min_cluster_index = -1; float min_cluster_distance = 99999; bool flip = false; for (unsigned int k=0; k v = C.at(k).h / C.at(k).n; bool f = false; float d = 0; for (auto m : m_Metrics) d += m->CalculateDistance(t, v, f); d /= m_Metrics.size(); if (d=0 && min_cluster_distance outC; //#pragma omp parallel for for (int c=0; c<(int)C.size(); c++) { std::vector< Cluster > tempC = ClusterStep(C.at(c).I, distances); AppendCluster(outC, tempC); } return outC; } else return C; } void TractClusteringFilter::AppendCluster(std::vector< Cluster >& a, std::vector< Cluster >&b) { for (auto c : b) a.push_back(c); } void TractClusteringFilter::MergeDuplicateClusters(std::vector< TractClusteringFilter::Cluster >& clusters) { if (m_MergeDuplicateThreshold<0) m_MergeDuplicateThreshold = m_Distances.at(0)/2; bool found = true; MITK_INFO << "Merging duplicate clusters with distance threshold " << m_MergeDuplicateThreshold; int start = 0; while (found && m_MergeDuplicateThreshold>mitk::eps) { std::cout << " \r"; std::cout << "Number of clusters: " << clusters.size() << '\r'; cout.flush(); found = false; for (int k1=start; k1<(int)clusters.size(); ++k1) { Cluster c1 = clusters.at(k1); vnl_matrix t = c1.h / c1.n; std::vector< int > merge_indices; std::vector< bool > flip_indices; #pragma omp parallel for for (int k2=0; k2<(int)clusters.size(); ++k2) { if (k1!=k2) { Cluster c2 = clusters.at(k2); vnl_matrix v = c2.h / c2.n; bool f = false; float d = 0; for (auto m : m_Metrics) d += m->CalculateDistance(t, v, f); d /= m_Metrics.size(); #pragma omp critical if (d TractClusteringFilter::AddToKnownClusters(std::vector< long > f_indices, std::vector >& centroids) { float dist_thres = m_Distances.at(0); int N = f_indices.size(); std::vector< Cluster > C; vnl_matrix zero_h; zero_h.set_size(T.at(0).rows(), T.at(0).cols()); zero_h.fill(0.0); Cluster no_fit; no_fit.h = zero_h; for (unsigned int i=0; i t = T.at(f_indices.at(i)); int min_cluster_index = -1; float min_cluster_distance = 99999; bool flip = false; if (CalcOverlap(t)>=m_OverlapThreshold) { int c_idx = 0; for (vnl_matrix centroid : centroids) { bool f = false; float d = 0; for (auto m : m_Metrics) d += m->CalculateDistance(t, centroid, f); d /= m_Metrics.size(); if (d=0 && min_cluster_distance f_indices; for (unsigned int i=0; i clusters; if (m_InCentroids.IsNull()) { MITK_INFO << "Clustering fibers"; clusters = ClusterStep(f_indices, m_Distances); MITK_INFO << "Number of clusters: " << clusters.size(); MergeDuplicateClusters(clusters); std::sort(clusters.begin(),clusters.end()); } else { std::vector > centroids = ResampleFibers(m_InCentroids); if (centroids.empty()) { MITK_INFO << "No fibers in centroid tractogram!"; return; } MITK_INFO << "Clustering with input centroids"; clusters = AddToKnownClusters(f_indices, centroids); no_match = clusters.back(); clusters.pop_back(); MITK_INFO << "Number of clusters: " << clusters.size(); MergeDuplicateClusters(clusters); } MITK_INFO << "Clustering finished"; int max = clusters.size()-1; if (m_MaxClusters>0 && clusters.size()-1>m_MaxClusters) max = m_MaxClusters; + int skipped = 0; for (int i=clusters.size()-1; i>=0; --i) { Cluster c = clusters.at(i); if ( c.n>=(int)m_MinClusterSize && !(m_MaxClusters>0 && clusters.size()-i>m_MaxClusters) ) { m_OutClusters.push_back(c); vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = m_Tractogram->GeneratePolyDataByIds(c.I, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); if (max>0) fib->SetFiberWeights((float)i/max); m_OutTractograms.push_back(fib); m_OutFiberIndices.push_back(c.I); // create centroid vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer polyData = vtkSmartPointer::New(); vtkSmartPointer container = vtkSmartPointer::New(); vnl_matrix centroid_points = c.h / c.n; for (unsigned int j=0; jInsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); polyData->SetPoints(vtkNewPoints); polyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer centroid = mitk::FiberBundle::New(polyData); centroid->SetFiberColors(255, 255, 255); m_OutCentroids.push_back(centroid); } + else + { + skipped++; + } } - MITK_INFO << "Final number of clusters: " << m_OutTractograms.size(); + MITK_INFO << "Final number of clusters: " << m_OutTractograms.size() << " (discarded " << skipped << " clusters)"; int w = 0; for (auto fib : m_OutTractograms) { if (m_OutTractograms.size()>1) { fib->SetFiberWeights((float)w/(m_OutTractograms.size()-1)); m_OutCentroids.at(w)->SetFiberWeights((float)w/(m_OutTractograms.size()-1)); } else { fib->SetFiberWeights(1); m_OutCentroids.at(w)->SetFiberWeights(1); } fib->ColorFibersByFiberWeights(false, false); ++w; } if (no_match.n>0) { vtkSmartPointer weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = m_Tractogram->GeneratePolyDataByIds(no_match.I, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); fib->SetFiberColors(0, 0, 0); m_OutFiberIndices.push_back(no_match.I); m_OutTractograms.push_back(fib); } } } diff --git a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp index 2c8ac8ca6c..3b77f7fa17 100755 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp @@ -1,2546 +1,2544 @@ /*=================================================================== 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 #include const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs"; 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; else { this->m_FiberPolyData->SetPoints(vtkSmartPointer::New()); this->m_FiberPolyData->SetLines(vtkSmartPointer::New()); } this->UpdateFiberGeometry(); this->GenerateFiberIds(); this->ColorFibersByOrientation(); } mitk::FiberBundle::~FiberBundle() { } mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy() { mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); newFib->SetFiberColors(this->m_FiberColors); newFib->SetFiberWeights(this->m_FiberWeights); return newFib; } vtkSmartPointer mitk::FiberBundle::GeneratePolyDataByIds(std::vector fiberIds, vtkSmartPointer weights) { vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); weights->SetNumberOfValues(fiberIds.size()); int counter = 0; auto finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*finIt < 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]); } weights->InsertValue(counter, this->GetFiberWeight(*finIt)); newLineSet->InsertNextCell(newFiber); ++finIt; ++counter; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundles(std::vector< mitk::FiberBundle::Pointer > fibs) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); int num_weights = this->GetNumFibers(); for (auto fib : fibs) num_weights += fib->GetNumFibers(); weights->SetNumberOfValues(num_weights); 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++; } for (auto fib : fibs) { // 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; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib) { if (fib==nullptr) return this->GetDeepCopy(); MITK_INFO << "Adding fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); unsigned int counter = 0; for (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; } // Only retain fibers with a weight larger than the specified threshold mitk::FiberBundle::Pointer mitk::FiberBundle::FilterByWeights(float weight_thr, bool invert) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector weights; for (int i=0; iGetNumFibers(); i++) { if ( (invert && this->GetFiberWeight(i)>weight_thr) || (!invert && this->GetFiberWeight(i)<=weight_thr)) continue; vtkCell* cell = m_FiberPolyData->GetCell(i); 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); } vNewLines->InsertNextCell(container); weights.push_back(this->GetFiberWeight(i)); } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); for (unsigned int i=0; iSetFiberWeight(i, weights.at(i)); return newFib; } // Only retain a subsample of the fibers mitk::FiberBundle::Pointer mitk::FiberBundle::SubsampleFibers(float factor) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); int new_num_fibs = this->GetNumFibers()*factor; MITK_INFO << "Subsampling fibers with factor " << factor << "(" << new_num_fibs << "/" << this->GetNumFibers() << ")"; // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(new_num_fibs); std::vector< int > ids; for (int i=0; iGetNumFibers(); i++) ids.push_back(i); std::random_shuffle(ids.begin(), ids.end()); unsigned int counter = 0; for (int i=0; iGetCell(ids.at(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(ids.at(i))); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // subtract two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib) { if (fib==nullptr) return this->GetDeepCopy(); MITK_INFO << "Subtracting fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector< std::vector< itk::Point > > points1; for( 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 mitk::FiberBundle::New(); // 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); std::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; + 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; + 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++; } } if (inside+outside==0) outside = 1; return (float)inside/(inside+outside); } 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 weights = vtkSmartPointer::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp, weights); mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(pTmp); fib->SetFiberWeights(weights); return fib; } std::vector mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage) { std::vector result; if (roi==nullptr || roi->GetData()==nullptr) return result; mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); if (!pfc.IsNull()) // handle composite { DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); if (children->size()==0) return result; switch (pfc->getOperationType()) { case 0: // AND { MITK_INFO << "AND"; result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(std::min(result.size(),inRoi.size())); it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( 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(itk::ScalableAffineTransform< mitk::ScalarType >::Pointer transform) { 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; j p = GetItkPoint(points->GetPoint(j)); p = transform->TransformPoint(p); vtkIdType id = vtkNewPoints->InsertNextPoint(p.GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) { 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 MITK_INFO << "Smoothing fibers"; vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); std::vector< vtkSmartPointer > resampled_streamlines; resampled_streamlines.resize(m_NumFibers); boost::progress_display disp(m_NumFibers); #pragma omp parallel for for (int i=0; i newPoints = vtkSmartPointer::New(); float length = 0; #pragma omp critical { length = m_FiberLengths.at(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(); #pragma omp critical { for (int j=0; jGetNumberOfPoints(); j++) { vtkIdType id = vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); smoothLine->GetPointIds()->InsertNextId(id); } resampled_streamlines[i] = smoothLine; } } for (auto container : resampled_streamlines) { vtkSmoothCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ResampleSpline(float pointDistance) { ResampleSpline(pointDistance, 0, 0, 0 ); } unsigned 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::ResampleToNumPoints(unsigned int targetPoints) { if (targetPoints<2) mitkThrow() << "Minimum two points required for resampling!"; MITK_INFO << "Resampling fibers (number of points " << targetPoints << ")"; bool unequal_fibs = true; while (unequal_fibs) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); unequal_fibs = false; //#pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; double seg_len = 0; //#pragma omp critical { weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); if ((unsigned int)numPoints!=targetPoints) seg_len = this->GetFiberLength(i)/(targetPoints-1);; vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); //#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 >= seg_len && seg_len>0) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-seg_len <= mitk::eps ) { vec.normalize(); newV += vec * seg_len; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - seg_len*seg_len; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } //#pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if ( (j==vertices.size()-1 && new_dist>0.0001) || seg_len==0) { //#pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } //#pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); vtkNewCells->InsertNextCell(container); if (container->GetNumberOfPoints()!=targetPoints) unequal_fibs = true; } } if (vtkNewCells->GetNumberOfCells()>0) { SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } } void mitk::FiberBundle::ResampleLinear(double pointDistance) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Resampling fibers (linear)"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); std::vector< vtkSmartPointer > resampled_streamlines; resampled_streamlines.resize(m_FiberPolyData->GetNumberOfCells()); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; #pragma omp critical { ++disp; 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 { resampled_streamlines[i] = container; } } for (auto container : resampled_streamlines) { vtkNewCells->InsertNextCell(container); } if (vtkNewCells->GetNumberOfCells()>0) { m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } // reapply selected colorcoding in case PolyData structure has changed bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps) { if (fib==nullptr) { MITK_INFO << "Reference bundle is nullptr!"; return false; } if (m_NumFibers!=fib->GetNumFibers()) { MITK_INFO << "Unequal number of fibers!"; MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); return false; } for (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/Testing/mitkFiberProcessingTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberProcessingTest.cpp index a1659473f4..0b245716c1 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberProcessingTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkFiberProcessingTest.cpp @@ -1,302 +1,304 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include #include #include #include #include #include #include #include "mitkTestFixture.h" class mitkFiberProcessingTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkFiberProcessingTestSuite); MITK_TEST(Test1); MITK_TEST(Test2); MITK_TEST(Test3); MITK_TEST(Test4); MITK_TEST(Test5); MITK_TEST(Test6); MITK_TEST(Test7); MITK_TEST(Test8); MITK_TEST(Test9); MITK_TEST(Test10); MITK_TEST(Test11); MITK_TEST(Test12); MITK_TEST(Test13); MITK_TEST(Test14); MITK_TEST(Test15); MITK_TEST(Test16); MITK_TEST(Test17); MITK_TEST(Test18); CPPUNIT_TEST_SUITE_END(); typedef itk::Image ItkUcharImgType; private: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ mitk::FiberBundle::Pointer original; ItkUcharImgType::Pointer mask; public: void setUp() override { omp_set_num_threads(1); original = nullptr; original = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/original.fib")).front().GetPointer()); mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/MASK.nrrd"))[0].GetPointer()); mask = ItkUcharImgType::New(); mitk::CastToItkImage(img, mask); } void tearDown() override { original = nullptr; } void Test1() { MITK_INFO << "TEST 1: Remove by direction"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); vnl_vector_fixed dir; dir[0] = 0; dir[1] = 1; dir[2] = 0; fib->RemoveDir(dir,cos(5.0*M_PI/180)); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/remove_direction.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test2() { MITK_INFO << "TEST 2: Remove by length"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->RemoveShortFibers(30); fib->RemoveLongFibers(40); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/remove_length.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test3() { MITK_INFO << "TEST 3: Remove by curvature 1"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); auto filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(30); filter->SetDistance(5); filter->SetRemoveFibers(false); + filter->SetUseMedian(true); filter->Update(); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/remove_curvature1.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(filter->GetOutputFiberBundle())); } void Test4() { MITK_INFO << "TEST 4: Remove by curvature 2"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); auto filter = itk::FiberCurvatureFilter::New(); filter->SetInputFiberBundle(fib); filter->SetAngularDeviation(30); filter->SetDistance(5); filter->SetRemoveFibers(true); + filter->SetUseMedian(true); filter->Update(); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/remove_curvature2.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(filter->GetOutputFiberBundle())); } void Test5() { MITK_INFO << "TEST 5: Remove outside mask"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib = fib->RemoveFibersOutside(mask); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/remove_outside.fib")).front().GetPointer()); mitk::IOUtil::Save(fib, mitk::IOUtil::GetTempPath()+"remove_outside.fib"); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test6() { MITK_INFO << "TEST 6: Remove inside mask"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib = fib->RemoveFibersOutside(mask, true); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/remove_inside.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test7() { MITK_INFO << "TEST 7: Modify resample spline"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->ResampleSpline(5); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/modify_resample.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test8() { MITK_INFO << "TEST 8: Modify compress"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->Compress(0.1); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/modify_compress.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test9() { MITK_INFO << "TEST 9: Modify sagittal mirror"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->MirrorFibers(0); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/modify_sagittal_mirror.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test10() { MITK_INFO << "TEST 10: Modify coronal mirror"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->MirrorFibers(1); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/modify_coronal_mirror.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test11() { MITK_INFO << "TEST 11: Modify axial mirror"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->MirrorFibers(2); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/modify_axial_mirror.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test12() { MITK_INFO << "TEST 12: Weight and join"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->SetFiberWeights(0.1); mitk::FiberBundle::Pointer fib2 = original->GetDeepCopy(); fib2->SetFiberWeight(3, 0.5); fib = fib->AddBundle(fib2); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/modify_weighted_joined.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Number of fibers", ref->GetNumFibers() == fib->GetNumFibers()); for (int i=0; iGetNumFibers(); i++) CPPUNIT_ASSERT_MESSAGE("Fiber weights", ref->GetFiberWeight(i) == fib->GetFiberWeight(i)); } void Test13() { MITK_INFO << "TEST 13: Subtract"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); mitk::FiberBundle::Pointer fib2 = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/remove_length.fib")).front().GetPointer()); fib = fib->SubtractBundle(fib2); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/subtracted.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test14() { MITK_INFO << "TEST 14: rotate"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->TransformFibers(1,2,3,0,0,0); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/transform_rotate.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test15() { MITK_INFO << "TEST 15: translate"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->TransformFibers(0,0,0,1,2,3); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/transform_translate.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test16() { MITK_INFO << "TEST 16: scale 1"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->ScaleFibers(0.1, 0.2, 0.3); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/transform_scale1.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test17() { MITK_INFO << "TEST 17: scale 2"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->ScaleFibers(0.1, 0.2, 0.3, false); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/transform_scale2.fib")).front().GetPointer()); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } void Test18() { MITK_INFO << "TEST 18: Modify resample linear"; mitk::FiberBundle::Pointer fib = original->GetDeepCopy(); fib->ResampleLinear(); mitk::FiberBundle::Pointer ref = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/FiberProcessing/modify_resample_linear.fib")).front().GetPointer()); mitk::IOUtil::Save(fib, mitk::IOUtil::GetTempPath()+"modify_resample_linear.fib"); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(fib)); } }; MITK_TEST_SUITE_REGISTRATION(mitkFiberProcessing) diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp index 36f231a281..b83b593cad 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp @@ -1,109 +1,108 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include #include #include #include #include #include #include #include #include #include #include #include "mitkTestFixture.h" class mitkMachineLearningTrackingTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkMachineLearningTrackingTestSuite); MITK_TEST(Track1); CPPUNIT_TEST_SUITE_END(); typedef itk::Image ItkFloatImgType; private: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ mitk::FiberBundle::Pointer ref; mitk::TrackingHandlerRandomForest<6, 100>* tfh; mitk::Image::Pointer dwi; ItkFloatImgType::Pointer seed; public: void setUp() override { ref = nullptr; tfh = new mitk::TrackingHandlerRandomForest<6,100>(); std::vector fibInfile = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/ReferenceTracts.fib")); mitk::BaseData::Pointer baseData = fibInfile.at(0); ref = dynamic_cast(baseData.GetPointer()); dwi = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/DiffusionImage.dwi"))[0].GetPointer()); mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/seed.nrrd"))[0].GetPointer()); seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); mitk::TractographyForest::Pointer forest = dynamic_cast(mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/MachineLearningTracking/forest.rf"))[0].GetPointer()); tfh->SetForest(forest); tfh->AddDwi(dwi); } void tearDown() override { delete tfh; ref = nullptr; } void Track1() { omp_set_num_threads(1); typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); tracker->SetDemoMode(false); tracker->SetInterpolateMask(false); tracker->SetSeedImage(seed); tracker->SetSeedsPerVoxel(1); tracker->SetStepSize(-1); tracker->SetAngularThreshold(45); tracker->SetMinTractLength(20); tracker->SetMaxTractLength(400); tracker->SetTrackingHandler(tfh); - tracker->SetAposterioriCurvCheck(false); tracker->SetAvoidStop(true); tracker->SetSamplingDistance(0.5); tracker->SetRandomSampling(false); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); //MITK_INFO << mitk::IOUtil::GetTempPath() << "ReferenceTracts.fib"; //mitk::IOUtil::Save(outFib, mitk::IOUtil::GetTempPath()+"ReferenceTracts.fib"); CPPUNIT_ASSERT_MESSAGE("Should be equal", ref->Equals(outFib)); } }; MITK_TEST_SUITE_REGISTRATION(mitkMachineLearningTracking) diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp index 9e23f45457..d2db4e6d9e 100755 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp @@ -1,427 +1,426 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include class mitkStreamlineTractographyTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkStreamlineTractographyTestSuite); MITK_TEST(Test_Peak1); MITK_TEST(Test_Peak2); MITK_TEST(Test_Tensor1); MITK_TEST(Test_Tensor2); MITK_TEST(Test_Tensor3); MITK_TEST(Test_Odf1); MITK_TEST(Test_Odf2); MITK_TEST(Test_Odf3); MITK_TEST(Test_Odf4); MITK_TEST(Test_Odf5); MITK_TEST(Test_Odf6); CPPUNIT_TEST_SUITE_END(); typedef itk::VectorImage< short, 3> ItkDwiType; private: public: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ typedef itk::Image ItkFloatImgType; mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itk_odf_image; mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itk_tensor_image; mitk::TrackingHandlerPeaks::PeakImgType::Pointer itk_peak_image; ItkFloatImgType::Pointer itk_seed_image; ItkFloatImgType::Pointer itk_mask_image; ItkFloatImgType::Pointer itk_gfa_image; float gfa_threshold; float odf_threshold; float peak_threshold; itk::StreamlineTrackingFilter::Pointer tracker; void setUp() override { omp_set_num_threads(1); gfa_threshold = 0.2; odf_threshold = 0.1; peak_threshold = 0.1; mitk::Image::Pointer odf_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/qball_image.qbi")); mitk::Image::Pointer tensor_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/tensor_image.dti")); mitk::Image::Pointer peak_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/qball_peak_image.nii.gz")); mitk::Image::Pointer seed_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/seed_image.nii.gz")); mitk::Image::Pointer mask_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/mask_image.nii.gz")); mitk::Image::Pointer gfa_image = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/gfa_image.nii.gz")); { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(peak_image); caster->Update(); itk_peak_image = caster->GetOutput(); } { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(tensor_image); caster->Update(); itk_tensor_image = caster->GetOutput(); } { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(odf_image); caster->Update(); itk_odf_image = caster->GetOutput(); } itk_gfa_image = ItkFloatImgType::New(); mitk::CastToItkImage(gfa_image, itk_gfa_image); itk_seed_image = ItkFloatImgType::New(); mitk::CastToItkImage(seed_image, itk_seed_image); itk_mask_image = ItkFloatImgType::New(); mitk::CastToItkImage(mask_image, itk_mask_image); } mitk::FiberBundle::Pointer LoadReferenceFib(std::string filename) { mitk::FiberBundle::Pointer fib = nullptr; if (itksys::SystemTools::FileExists(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename))) { mitk::BaseData::Pointer baseData = mitk::IOUtil::Load(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename)).at(0); fib = dynamic_cast(baseData.GetPointer()); } return fib; } mitk::Image::Pointer LoadReferenceImage(std::string filename) { mitk::Image::Pointer img = nullptr; if (itksys::SystemTools::FileExists(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename))) { img = mitk::IOUtil::LoadImage(GetTestDataFilePath("DiffusionImaging/StreamlineTractography/ReferenceFibs/" + filename)); } return img; } void SetupTracker(mitk::TrackingDataHandler* handler) { tracker = itk::StreamlineTrackingFilter::New(); tracker->SetRandom(false); tracker->SetInterpolateMask(false); tracker->SetNumberOfSamples(0); tracker->SetAngularThreshold(-1); tracker->SetMaskImage(itk_mask_image); tracker->SetSeedImage(itk_seed_image); tracker->SetStoppingRegions(nullptr); tracker->SetSeedsPerVoxel(1); tracker->SetStepSize(0.5); tracker->SetSamplingDistance(0.25); tracker->SetUseStopVotes(true); tracker->SetOnlyForwardSamples(true); - tracker->SetAposterioriCurvCheck(false); tracker->SetMinTractLength(20); tracker->SetMaxNumTracts(-1); tracker->SetTrackingHandler(handler); tracker->SetUseOutputProbabilityMap(false); } void tearDown() override { } void CheckFibResult(std::string ref_file, mitk::FiberBundle::Pointer test_fib) { mitk::FiberBundle::Pointer ref = LoadReferenceFib(ref_file); if (ref.IsNull()) { mitk::IOUtil::Save(test_fib, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Reference file not found. Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } else { bool is_equal = ref->Equals(test_fib); if (!is_equal) { mitk::IOUtil::Save(test_fib, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Tractograms are not equal! Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } } } void CheckImageResult(std::string ref_file, mitk::Image::Pointer test_img) { mitk::Image::Pointer ref = LoadReferenceImage(ref_file); if (ref.IsNull()) { mitk::IOUtil::Save(test_img, mitk::IOUtil::GetTempPath()+ref_file); CPPUNIT_FAIL("Reference file not found. Saving test file to " + mitk::IOUtil::GetTempPath() + ref_file); } else { MITK_ASSERT_EQUAL(test_img, ref, "Images should be equal"); } } void Test_Peak1() { mitk::TrackingHandlerPeaks* handler = new mitk::TrackingHandlerPeaks(); handler->SetPeakImage(itk_peak_image); handler->SetPeakThreshold(peak_threshold); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Peak1.fib", outFib); delete handler; } void Test_Peak2() { mitk::TrackingHandlerPeaks* handler = new mitk::TrackingHandlerPeaks(); handler->SetPeakImage(itk_peak_image); handler->SetPeakThreshold(peak_threshold); handler->SetInterpolate(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Peak2.fib", outFib); delete handler; } void Test_Tensor1() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); handler->SetFaThreshold(gfa_threshold); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor1.fib", outFib); delete handler; } void Test_Tensor2() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); handler->SetFaThreshold(gfa_threshold); handler->SetInterpolate(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor2.fib", outFib); delete handler; } void Test_Tensor3() { mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetTensorImage(itk_tensor_image); handler->SetFaThreshold(gfa_threshold); handler->SetInterpolate(false); handler->SetF(0); handler->SetG(1); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Tensor3.fib", outFib); delete handler; } void Test_Odf1() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf1.fib", outFib); delete handler; } void Test_Odf2() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf2.fib", outFib); delete handler; } void Test_Odf3() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); handler->SetInterpolate(false); SetupTracker(handler); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf3.fib", outFib); delete handler; } void Test_Odf4() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); SetupTracker(handler); tracker->SetSeedsPerVoxel(3); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf4.fib", outFib); delete handler; } void Test_Odf5() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); SetupTracker(handler); tracker->SetSeedsPerVoxel(3); tracker->Update(); vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); CheckFibResult("Test_Odf5.fib", outFib); delete handler; } void Test_Odf6() { mitk::TrackingHandlerOdf* handler = new mitk::TrackingHandlerOdf(); handler->SetOdfImage(itk_odf_image); handler->SetGfaThreshold(gfa_threshold); handler->SetOdfThreshold(0); handler->SetSharpenOdfs(true); handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); SetupTracker(handler); tracker->SetSeedsPerVoxel(10); tracker->SetUseOutputProbabilityMap(true); tracker->Update(); itk::StreamlineTrackingFilter::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::IOUtil::Save(img, mitk::IOUtil::GetTempPath()+"Test_Odf6.nrrd"); CheckImageResult("Test_Odf6.nrrd", img); delete handler; } }; MITK_TEST_SUITE_REGISTRATION(mitkStreamlineTractography) diff --git a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp index bee4bedf2e..39f0b36536 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/FiberClustering.cpp @@ -1,228 +1,233 @@ /*=================================================================== 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 mitk::FiberBundle::Pointer LoadFib(std::string filename) { std::vector fibInfile = mitk::IOUtil::Load(filename); if( fibInfile.empty() ) std::cout << "File " << filename << " could not be read!"; mitk::BaseData::Pointer baseData = fibInfile.at(0); return dynamic_cast(baseData.GetPointer()); } /*! \brief Spatially cluster fibers */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Fiber Clustering"); parser.setCategory("Fiber Processing"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("", "i", mitkCommandLineParser::InputFile, "Input:", "input fiber bundle (.fib, .trk, .tck)", us::Any(), false); parser.addArgument("", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output root", us::Any(), false); parser.addArgument("cluster_size", "", mitkCommandLineParser::Int, "Cluster size:", "", 10); parser.addArgument("fiber_points", "", mitkCommandLineParser::Int, "Fiber points:", "", 12); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. fibers per cluster:", "", 1); parser.addArgument("max_clusters", "", mitkCommandLineParser::Int, "Max. clusters:", ""); parser.addArgument("merge_clusters", "", mitkCommandLineParser::Float, "Merge clusters:", "", -1.0); parser.addArgument("output_centroids", "", mitkCommandLineParser::Bool, "Output centroids:", ""); parser.addArgument("metrics", "", mitkCommandLineParser::StringList, "Metrics:", "EU_MEAN, EU_STD, EU_MAX, ANAT, MAP, ANGLES"); parser.addArgument("input_centroids", "", mitkCommandLineParser::String, "Input centroids:", ""); parser.addArgument("scalar_map", "", mitkCommandLineParser::String, "Scalar map:", ""); parser.addArgument("parcellation", "", mitkCommandLineParser::String, "Parcellation:", ""); + parser.addArgument("file_ending", "", mitkCommandLineParser::String, "File ending:", ""); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFileName = us::any_cast(parsedArgs["i"]); std::string out_root = us::any_cast(parsedArgs["o"]); int cluster_size = 10; if (parsedArgs.count("cluster_size")) cluster_size = us::any_cast(parsedArgs["cluster_size"]); int fiber_points = 12; if (parsedArgs.count("fiber_points")) fiber_points = us::any_cast(parsedArgs["fiber_points"]); int min_fibers = 1; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); int max_clusters = 0; if (parsedArgs.count("max_clusters")) max_clusters = us::any_cast(parsedArgs["max_clusters"]); float merge_clusters = -1.0; if (parsedArgs.count("merge_clusters")) merge_clusters = us::any_cast(parsedArgs["merge_clusters"]); bool output_centroids = false; if (parsedArgs.count("output_centroids")) output_centroids = us::any_cast(parsedArgs["output_centroids"]); std::vector< std::string > metric_strings = {"EU_MEAN"}; if (parsedArgs.count("metrics")) metric_strings = us::any_cast(parsedArgs["metrics"]); std::string input_centroids = ""; if (parsedArgs.count("input_centroids")) input_centroids = us::any_cast(parsedArgs["input_centroids"]); std::string scalar_map = ""; if (parsedArgs.count("scalar_map")) scalar_map = us::any_cast(parsedArgs["scalar_map"]); std::string parcellation = ""; if (parsedArgs.count("parcellation")) parcellation = us::any_cast(parsedArgs["parcellation"]); + std::string file_ending = ".fib"; + if (parsedArgs.count("file_ending")) + file_ending = us::any_cast(parsedArgs["file_ending"]); + try { typedef itk::Image< float, 3 > FloatImageType; typedef itk::Image< short, 3 > ShortImageType; mitk::FiberBundle::Pointer fib = LoadFib(inFileName); float max_d = 0; int i=1; std::vector< float > distances; while (max_d < fib->GetGeometry()->GetDiagonalLength()/2) { distances.push_back(cluster_size*i); max_d = cluster_size*i; ++i; } itk::TractClusteringFilter::Pointer clusterer = itk::TractClusteringFilter::New(); clusterer->SetDistances(distances); clusterer->SetTractogram(fib); if (input_centroids!="") { mitk::FiberBundle::Pointer in_centroids = LoadFib(input_centroids); clusterer->SetInCentroids(in_centroids); } std::vector< mitk::ClusteringMetric* > metrics; for (auto m : metric_strings) { MITK_INFO << "Metric: " << m; if (m=="EU_MEAN") metrics.push_back({new mitk::ClusteringMetricEuclideanMean()}); else if (m=="EU_STD") metrics.push_back({new mitk::ClusteringMetricEuclideanStd()}); else if (m=="EU_MAX") metrics.push_back({new mitk::ClusteringMetricEuclideanMax()}); else if (m=="ANGLES") metrics.push_back({new mitk::ClusteringMetricInnerAngles()}); else if (m=="MAP" && scalar_map!="") { mitk::Image::Pointer mitk_map = dynamic_cast(mitk::IOUtil::Load(scalar_map)[0].GetPointer()); if (mitk_map->GetDimension()==3) { FloatImageType::Pointer itk_map = FloatImageType::New(); mitk::CastToItkImage(mitk_map, itk_map); mitk::ClusteringMetricScalarMap* metric = new mitk::ClusteringMetricScalarMap(); metric->SetImages({itk_map}); metric->SetScale(distances.at(0)); metrics.push_back(metric); } } else if (m=="ANAT" && parcellation!="") { mitk::Image::Pointer mitk_map = dynamic_cast(mitk::IOUtil::Load(parcellation)[0].GetPointer()); if (mitk_map->GetDimension()==3) { ShortImageType::Pointer itk_map = ShortImageType::New(); mitk::CastToItkImage(mitk_map, itk_map); mitk::ClusteringMetricAnatomic* metric = new mitk::ClusteringMetricAnatomic(); metric->SetParcellations({itk_map}); metrics.push_back(metric); } } } if (metrics.empty()) { MITK_INFO << "No metric selected!"; return EXIT_FAILURE; } clusterer->SetMetrics(metrics); clusterer->SetMergeDuplicateThreshold(merge_clusters); clusterer->SetNumPoints(fiber_points); clusterer->SetMaxClusters(max_clusters); clusterer->SetMinClusterSize(min_fibers); clusterer->Update(); std::vector tracts = clusterer->GetOutTractograms(); std::vector centroids = clusterer->GetOutCentroids(); MITK_INFO << "Saving clusters"; std::streambuf *old = cout.rdbuf(); // <-- save std::stringstream ss; std::cout.rdbuf (ss.rdbuf()); // <-- redirect unsigned int c = 0; for (auto f : tracts) { - mitk::IOUtil::Save(f, out_root + "Cluster_" + boost::lexical_cast(c) + ".fib"); + mitk::IOUtil::Save(f, out_root + "Cluster_" + boost::lexical_cast(c) + file_ending); if (output_centroids) - mitk::IOUtil::Save(centroids.at(c), out_root + "Centroid_" + boost::lexical_cast(c) + ".fib"); + mitk::IOUtil::Save(centroids.at(c), out_root + "Centroid_" + boost::lexical_cast(c) + file_ending); ++c; } std::cout.rdbuf (old); // <-- restore } 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/FiberProcessing/TractDensity.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensity.cpp index a34fd5a0ac..0c2e861845 100644 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensity.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensity.cpp @@ -1,201 +1,206 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include mitk::FiberBundle::Pointer LoadFib(std::string filename) { - std::vector fibInfile = mitk::IOUtil::Load(filename); - if( fibInfile.empty() ) - std::cout << "File " << filename << " could not be read!"; - mitk::BaseData::Pointer baseData = fibInfile.at(0); - return dynamic_cast(baseData.GetPointer()); + std::vector fibInfile = mitk::IOUtil::Load(filename); + if( fibInfile.empty() ) + std::cout << "File " << filename << " could not be read!"; + mitk::BaseData::Pointer baseData = fibInfile.at(0); + return dynamic_cast(baseData.GetPointer()); } /*! \brief Modify input tractogram: fiber resampling, compression, pruning and transformation. */ int main(int argc, char* argv[]) { - mitkCommandLineParser parser; + mitkCommandLineParser parser; - parser.setTitle("Tract Density"); - parser.setCategory("Fiber Tracking and Processing Methods"); - parser.setDescription("Generate tract density image, fiber envelope or fiber endpoints image."); - parser.setContributor("MIC"); + parser.setTitle("Tract Density"); + parser.setCategory("Fiber Tracking and Processing Methods"); + parser.setDescription("Generate tract density image, fiber envelope or fiber endpoints image."); + parser.setContributor("MIC"); - parser.setArgumentPrefix("--", "-"); - parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle (.fib)", us::Any(), false); - parser.addArgument("output", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false); - parser.addArgument("binary", "", mitkCommandLineParser::Bool, "Binary output:", "calculate binary tract envelope", us::Any()); - parser.addArgument("endpoints", "", mitkCommandLineParser::Bool, "Output endpoints image:", "calculate image of fiber endpoints instead of mask", us::Any()); - parser.addArgument("reference_image", "", mitkCommandLineParser::String, "Reference image:", "output image will have geometry of this reference image", us::Any()); + parser.setArgumentPrefix("--", "-"); + parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input fiber bundle (.fib)", us::Any(), false); + parser.addArgument("output", "o", mitkCommandLineParser::String, "Output:", "output image", us::Any(), false); + parser.addArgument("binary", "", mitkCommandLineParser::Bool, "Binary output:", "calculate binary tract envelope", us::Any()); + parser.addArgument("normalize", "", mitkCommandLineParser::Bool, "Normalized output:", "normalize output to 0-1", us::Any()); + parser.addArgument("endpoints", "", mitkCommandLineParser::Bool, "Output endpoints image:", "calculate image of fiber endpoints instead of mask", us::Any()); + parser.addArgument("reference_image", "", mitkCommandLineParser::String, "Reference image:", "output image will have geometry of this reference image", us::Any()); - std::map parsedArgs = parser.parseArguments(argc, argv); - if (parsedArgs.size()==0) - return EXIT_FAILURE; + std::map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size()==0) + return EXIT_FAILURE; - bool binary = false; - if (parsedArgs.count("binary")) - binary = us::any_cast(parsedArgs["binary"]); + bool binary = false; + if (parsedArgs.count("binary")) + binary = us::any_cast(parsedArgs["binary"]); - bool endpoints = false; - if (parsedArgs.count("endpoints")) - endpoints = us::any_cast(parsedArgs["endpoints"]); + bool endpoints = false; + if (parsedArgs.count("endpoints")) + endpoints = us::any_cast(parsedArgs["endpoints"]); - std::string reference_image = ""; - if (parsedArgs.count("reference_image")) - reference_image = us::any_cast(parsedArgs["reference_image"]); + bool normalize = false; + if (parsedArgs.count("normalize")) + normalize = us::any_cast(parsedArgs["normalize"]); - std::string inFileName = us::any_cast(parsedArgs["input"]); - std::string outFileName = us::any_cast(parsedArgs["output"]); + std::string reference_image = ""; + if (parsedArgs.count("reference_image")) + reference_image = us::any_cast(parsedArgs["reference_image"]); - try - { - mitk::FiberBundle::Pointer fib = LoadFib(inFileName); - - mitk::Image::Pointer ref_img; - MITK_INFO << reference_image; - if (!reference_image.empty()) - ref_img = dynamic_cast(mitk::IOUtil::Load(reference_image)[0].GetPointer()); - - if (endpoints) - { - typedef unsigned int OutPixType; - typedef itk::Image OutImageType; - - typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; - ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); - generator->SetFiberBundle(fib); - - if (ref_img.IsNotNull()) - { - OutImageType::Pointer itkImage = OutImageType::New(); - CastToItkImage(ref_img, itkImage); - generator->SetInputImage(itkImage); - generator->SetUseImageGeometry(true); - - } - generator->Update(); - - // get output image - typedef itk::Image OutType; - OutType::Pointer outImg = generator->GetOutput(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - - mitk::IOUtil::Save(img, outFileName ); - } - else if (binary) - { - typedef unsigned char OutPixType; - typedef itk::Image OutImageType; - - itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); - generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); - generator->SetOutputAbsoluteValues(false); - generator->SetWorkOnFiberCopy(false); - - if (ref_img.IsNotNull()) - { - OutImageType::Pointer itkImage = OutImageType::New(); - CastToItkImage(ref_img, itkImage); - generator->SetInputImage(itkImage); - generator->SetUseImageGeometry(true); - - } - generator->Update(); - - // get output image - typedef itk::Image OutType; - OutType::Pointer outImg = generator->GetOutput(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - - mitk::IOUtil::Save(img, outFileName ); - } - else - { - typedef float OutPixType; - typedef itk::Image OutImageType; - - itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); - generator->SetFiberBundle(fib); - generator->SetBinaryOutput(binary); - generator->SetOutputAbsoluteValues(false); - generator->SetWorkOnFiberCopy(false); - - if (ref_img.IsNotNull()) - { - OutImageType::Pointer itkImage = OutImageType::New(); - CastToItkImage(ref_img, itkImage); - generator->SetInputImage(itkImage); - generator->SetUseImageGeometry(true); - - } - generator->Update(); - - // get output image - typedef itk::Image OutType; - OutType::Pointer outImg = generator->GetOutput(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - - mitk::IOUtil::Save(img, outFileName ); - } + std::string inFileName = us::any_cast(parsedArgs["input"]); + std::string outFileName = us::any_cast(parsedArgs["output"]); - } - catch (itk::ExceptionObject e) + try + { + mitk::FiberBundle::Pointer fib = LoadFib(inFileName); + + mitk::Image::Pointer ref_img; + MITK_INFO << reference_image; + if (!reference_image.empty()) + ref_img = dynamic_cast(mitk::IOUtil::Load(reference_image)[0].GetPointer()); + + if (endpoints) { - std::cout << e; - return EXIT_FAILURE; + typedef unsigned int OutPixType; + typedef itk::Image OutImageType; + + typedef itk::TractsToFiberEndingsImageFilter< OutImageType > ImageGeneratorType; + ImageGeneratorType::Pointer generator = ImageGeneratorType::New(); + generator->SetFiberBundle(fib); + + if (ref_img.IsNotNull()) + { + OutImageType::Pointer itkImage = OutImageType::New(); + CastToItkImage(ref_img, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + + } + generator->Update(); + + // get output image + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + mitk::IOUtil::Save(img, outFileName ); } - catch (std::exception e) + else if (binary) { - std::cout << e.what(); - return EXIT_FAILURE; + typedef unsigned char OutPixType; + typedef itk::Image OutImageType; + + itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); + generator->SetFiberBundle(fib); + generator->SetBinaryOutput(binary); + generator->SetOutputAbsoluteValues(!normalize); + generator->SetWorkOnFiberCopy(false); + + if (ref_img.IsNotNull()) + { + OutImageType::Pointer itkImage = OutImageType::New(); + CastToItkImage(ref_img, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + + } + generator->Update(); + + // get output image + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + mitk::IOUtil::Save(img, outFileName ); } - catch (...) + else { - std::cout << "ERROR!?!"; - return EXIT_FAILURE; + typedef float OutPixType; + typedef itk::Image OutImageType; + + itk::TractDensityImageFilter< OutImageType >::Pointer generator = itk::TractDensityImageFilter< OutImageType >::New(); + generator->SetFiberBundle(fib); + generator->SetBinaryOutput(binary); + generator->SetOutputAbsoluteValues(!normalize); + generator->SetWorkOnFiberCopy(false); + + if (ref_img.IsNotNull()) + { + OutImageType::Pointer itkImage = OutImageType::New(); + CastToItkImage(ref_img, itkImage); + generator->SetInputImage(itkImage); + generator->SetUseImageGeometry(true); + + } + generator->Update(); + + // get output image + typedef itk::Image OutType; + OutType::Pointer outImg = generator->GetOutput(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + mitk::IOUtil::Save(img, outFileName ); } - return EXIT_SUCCESS; + + } + 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/FiberProcessing/TractDensityFilter.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensityFilter.cpp index c99c28c7ef..d2b9a4b871 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensityFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/FiberProcessing/TractDensityFilter.cpp @@ -1,111 +1,111 @@ /*=================================================================== 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 "mitkCommandLineParser.h" #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include typedef itksys::SystemTools ist; typedef itk::Image ItkFloatImgType; /*! \brief Extract fibers from a tractogram using binary image ROIs */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Filter Outliers by Tract Density"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setContributor("MIC"); parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::String, "Input:", "input tractogram (.fib/.trk/.tck/.dcm)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::String, "Output:", "output tractogram", us::Any(), false); - parser.addArgument("threshold", "", mitkCommandLineParser::Float, "Threshold:", "positive means ROI image value threshold", false, 0.05); - parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap:", "positive means ROI image value threshold", false, 0.5); + parser.addArgument("threshold", "", mitkCommandLineParser::Float, "Threshold:", "positive means ROI image value threshold", 0.05, false); + parser.addArgument("overlap", "", mitkCommandLineParser::Float, "Overlap:", "positive means ROI image value threshold", 0.5, false); parser.addArgument("min_fibers", "", mitkCommandLineParser::Int, "Min. num. fibers:", "discard positive tracts with less fibers", 0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; std::string inFib = us::any_cast(parsedArgs["input"]); std::string outFib = us::any_cast(parsedArgs["out"]); int min_fibers = 0; if (parsedArgs.count("min_fibers")) min_fibers = us::any_cast(parsedArgs["min_fibers"]); float overlap = 0.5; if (parsedArgs.count("overlap")) overlap = us::any_cast(parsedArgs["overlap_fraction"]); float threshold = 0.05; if (parsedArgs.count("threshold")) threshold = us::any_cast(parsedArgs["threshold"]); try { mitk::FiberBundle::Pointer inputTractogram = dynamic_cast(mitk::IOUtil::Load(inFib)[0].GetPointer()); itk::TractDensityImageFilter< ItkFloatImgType >::Pointer generator = itk::TractDensityImageFilter< ItkFloatImgType >::New(); generator->SetFiberBundle(inputTractogram); generator->SetBinaryOutput(false); generator->SetOutputAbsoluteValues(false); generator->SetWorkOnFiberCopy(true); generator->Update(); itk::FiberExtractionFilter::Pointer extractor = itk::FiberExtractionFilter::New(); extractor->SetRoiImages({generator->GetOutput()}); extractor->SetInputFiberBundle(inputTractogram); extractor->SetOverlapFraction(overlap); extractor->SetInterpolate(true); extractor->SetThreshold(threshold); extractor->SetNoNegatives(true); extractor->Update(); if (extractor->GetPositives().at(0)->GetNumFibers()>=min_fibers) mitk::IOUtil::Save(extractor->GetPositives().at(0), outFib); } 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/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp index ce317deddf..bd84589eb5 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp @@ -1,456 +1,463 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false); parser.addArgument("algorithm", "a", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle/probability map", us::Any(), false); - parser.addArgument("stop_mask", "", mitkCommandLineParser::String, "Stop image:", "streamlines entering the binary mask will stop immediately", us::Any()); - parser.addArgument("target_image", "", mitkCommandLineParser::String, "Target image:", "streamlines not starting and ending in one of the regions in this image are discarded", us::Any()); + parser.addArgument("stop_image", "", mitkCommandLineParser::String, "Stop image:", "streamlines entering the binary mask will stop immediately", us::Any()); + parser.addArgument("target_image", "", mitkCommandLineParser::String, "Target image:", "a streamline is only considered valid if it starts and ends in the target and seed region region. without seed region, both endpoints need to be located in the target region", us::Any()); parser.addArgument("tracking_mask", "", mitkCommandLineParser::String, "Mask image:", "restrict tractography with a binary mask image", us::Any()); - parser.addArgument("seed_mask", "", mitkCommandLineParser::String, "Seed image:", "binary mask image defining seed voxels", us::Any()); + parser.addArgument("seed_image", "", mitkCommandLineParser::String, "Seed image:", "binary mask image defining seed voxels", us::Any()); parser.addArgument("sharpen_odfs", "", mitkCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); parser.addArgument("odf_cutoff", "", mitkCommandLineParser::Float, "ODF Cutoff:", "additional threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.1); parser.addArgument("step_size", "", mitkCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("angular_threshold", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size)"); parser.addArgument("min_tract_length", "", mitkCommandLineParser::Float, "Min. tract length:", "minimum fiber length (in mm)", 20); parser.addArgument("seeds", "", mitkCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); - parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded.", -1); + parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded", -1); + parser.addArgument("trials_per_seed", "", mitkCommandLineParser::Int, "Max. trials per seed:", "try each seed N times until a valid streamline is obtained (only for probabilistic)", 10); + parser.addArgument("loop_check", "", mitkCommandLineParser::Float, "Check for loops:", "threshold on angular stdev over the last 4 voxel lengths", -1); parser.addArgument("num_samples", "", mitkCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); - parser.addArgument("output_prob_map", "", mitkCommandLineParser::Bool, "Output probability map:", "output probability map instead of tractogram"); parser.addArgument("no_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate:", "don't interpolate image values"); parser.addArgument("flip_x", "", mitkCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); //parser.addArgument("apply_image_rotation", "", mitkCommandLineParser::Bool, "Apply image rotation:", "applies image rotation to image peaks (only for 'Peaks' algorithm). This is necessary in some cases, e.g. if the peaks were obtained with MRtrix."); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "Compress output fibers using the given error threshold (in mm)"); parser.addArgument("additional_images", "", mitkCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional Features)", us::Any()); // parameters for random forest based tractography parser.addArgument("forest", "", mitkCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any()); parser.addArgument("use_sh_features", "", mitkCommandLineParser::Bool, "Use SH features:", "use SH features"); // parameters for tensor tractography parser.addArgument("tend_f", "", mitkCommandLineParser::Float, "Weight f", "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkCommandLineParser::Float, "Weight g", "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["input"]); std::string outFile = us::any_cast(parsedArgs["out"]); std::string algorithm = us::any_cast(parsedArgs["algorithm"]); bool sharpen_odfs = false; if (parsedArgs.count("sharpen_odfs")) sharpen_odfs = us::any_cast(parsedArgs["sharpen_odfs"]); bool interpolate = true; if (parsedArgs.count("no_interpolation")) interpolate = !us::any_cast(parsedArgs["no_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); bool use_stop_votes = false; if (parsedArgs.count("use_stop_votes")) use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); bool use_only_forward_samples = false; if (parsedArgs.count("use_only_forward_samples")) use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); - bool output_prob_map = false; - if (parsedArgs.count("output_prob_map")) - output_prob_map = us::any_cast(parsedArgs["output_prob_map"]); - bool flip_x = false; if (parsedArgs.count("flip_x")) flip_x = us::any_cast(parsedArgs["flip_x"]); bool flip_y = false; if (parsedArgs.count("flip_y")) flip_y = us::any_cast(parsedArgs["flip_y"]); bool flip_z = false; if (parsedArgs.count("flip_z")) flip_z = us::any_cast(parsedArgs["flip_z"]); bool apply_image_rotation = false; if (parsedArgs.count("apply_image_rotation")) apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); float min_tract_length = 20; if (parsedArgs.count("min_tract_length")) min_tract_length = us::any_cast(parsedArgs["min_tract_length"]); + float loop_check = -1; + if (parsedArgs.count("loop_check")) + loop_check = us::any_cast(parsedArgs["loop_check"]); + std::string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); std::string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); std::string seedFile = ""; - if (parsedArgs.count("seed_mask")) - seedFile = us::any_cast(parsedArgs["seed_mask"]); + if (parsedArgs.count("seed_image")) + seedFile = us::any_cast(parsedArgs["seed_image"]); std::string targetFile = ""; if (parsedArgs.count("target_image")) targetFile = us::any_cast(parsedArgs["target_image"]); std::string stopFile = ""; - if (parsedArgs.count("stop_mask")) - stopFile = us::any_cast(parsedArgs["stop_mask"]); + if (parsedArgs.count("stop_image")) + stopFile = us::any_cast(parsedArgs["stop_image"]); float cutoff = 0.1; if (parsedArgs.count("cutoff")) cutoff = us::any_cast(parsedArgs["cutoff"]); float odf_cutoff = 0.1; if (parsedArgs.count("odf_cutoff")) odf_cutoff = us::any_cast(parsedArgs["odf_cutoff"]); float stepsize = -1; if (parsedArgs.count("step_size")) stepsize = us::any_cast(parsedArgs["step_size"]); float sampling_distance = -1; if (parsedArgs.count("sampling_distance")) sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); int num_samples = 0; if (parsedArgs.count("num_samples")) num_samples = us::any_cast(parsedArgs["num_samples"]); int seeds = 1; if (parsedArgs.count("seeds")) seeds = us::any_cast(parsedArgs["seeds"]); + unsigned int trials_per_seed = 10; + if (parsedArgs.count("trials_per_seed")) + trials_per_seed = us::any_cast(parsedArgs["trials_per_seed"]); + float tend_f = 1; if (parsedArgs.count("tend_f")) tend_f = us::any_cast(parsedArgs["tend_f"]); float tend_g = 0; if (parsedArgs.count("tend_g")) tend_g = us::any_cast(parsedArgs["tend_g"]); float angular_threshold = -1; if (parsedArgs.count("angular_threshold")) angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); unsigned int max_tracts = -1; if (parsedArgs.count("max_tracts")) max_tracts = us::any_cast(parsedArgs["max_tracts"]); - // LOAD DATASETS + std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); + if (ext != ".fib" && ext != ".trk") + { + MITK_INFO << "Output file format not supported. Use one of .fib, .trk, .nii, .nii.gz, .nrrd"; + return EXIT_FAILURE; + } + + // LOAD DATASETS mitkCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkFloatImgType; MITK_INFO << "loading input"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(mitk::IOUtil::Load(input_files.at(i))[0].GetPointer()); input_images.push_back(mitkImage); } ItkFloatImgType::Pointer mask; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(maskFile)[0].GetPointer()); mask = ItkFloatImgType::New(); mitk::CastToItkImage(img, mask); } ItkFloatImgType::Pointer seed; if (!seedFile.empty()) { MITK_INFO << "loading seed image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(seedFile)[0].GetPointer()); seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); } ItkFloatImgType::Pointer stop; if (!stopFile.empty()) { MITK_INFO << "loading stop image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(stopFile)[0].GetPointer()); stop = ItkFloatImgType::New(); mitk::CastToItkImage(img, stop); } ItkFloatImgType::Pointer target; if (!targetFile.empty()) { MITK_INFO << "loading target image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(targetFile)[0].GetPointer()); target = ItkFloatImgType::New(); mitk::CastToItkImage(img, target); } MITK_INFO << "loading additional images"; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(file)[0].GetPointer()); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); if (algorithm == "ProbTensor") { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkTensorImg = caster->GetOutput(); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkTensorImg ); filter->Update(); mitk::Image::Pointer image = mitk::Image::New(); FilterType::OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); input_images.clear(); input_images.push_back(image); sharpen_odfs = true; odf_cutoff = 0; } typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); mitk::TrackingDataHandler* handler; if (algorithm == "DetRF" || algorithm == "ProbRF") { mitk::TractographyForest::Pointer forest = dynamic_cast(mitk::IOUtil::Load(forestFile)[0].GetPointer()); if (forest.IsNull()) mitkThrow() << "Forest file " << forestFile << " could not be read."; if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->SetForest(forest); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } if (algorithm == "ProbRF") handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); } else if (algorithm == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetPeakImage(itkImg); dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); dynamic_cast(handler)->SetPeakThreshold(cutoff); } else if (algorithm == "DetTensor") { handler = new mitk::TrackingHandlerTensor(); for (auto input_image : input_images) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_image); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(handler)->AddTensorImage(itkImg); } dynamic_cast(handler)->SetFaThreshold(cutoff); dynamic_cast(handler)->SetF(tend_f); dynamic_cast(handler)->SetG(tend_g); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (algorithm == "DetODF" || algorithm == "ProbODF" || algorithm == "ProbTensor") { handler = new mitk::TrackingHandlerOdf(); typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetOdfImage(itkImg); dynamic_cast(handler)->SetGfaThreshold(cutoff); dynamic_cast(handler)->SetOdfThreshold(odf_cutoff); dynamic_cast(handler)->SetSharpenOdfs(sharpen_odfs); if (algorithm == "ProbODF" || algorithm == "ProbTensor") dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); if (algorithm == "ProbTensor") dynamic_cast(handler)->SetIsOdfFromTensor(true); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); } else { MITK_INFO << "Unknown tractography algorithm (" + algorithm+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } handler->SetInterpolate(interpolate); handler->SetFlipX(flip_x); handler->SetFlipY(flip_y); handler->SetFlipZ(flip_z); MITK_INFO << "Tractography algorithm: " << algorithm; - tracker->SetNumberOfSamples(num_samples); tracker->SetAngularThreshold(angular_threshold); tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetTargetRegions(target); tracker->SetSeedsPerVoxel(seeds); tracker->SetStepSize(stepsize); tracker->SetSamplingDistance(sampling_distance); tracker->SetUseStopVotes(use_stop_votes); tracker->SetOnlyForwardSamples(use_only_forward_samples); - tracker->SetAposterioriCurvCheck(false); + tracker->SetLoopCheck(loop_check); tracker->SetMaxNumTracts(max_tracts); + tracker->SetTrialsPerSeed(trials_per_seed); tracker->SetTrackingHandler(handler); - tracker->SetUseOutputProbabilityMap(output_prob_map); + if (ext != ".fib" && ext != ".trk") + tracker->SetUseOutputProbabilityMap(true); tracker->SetMinTractLength(min_tract_length); tracker->Update(); - std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); - - if (!output_prob_map) + if (ext == ".fib" || ext == ".trk") { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); - - if (ext != ".fib" && ext != ".trk" && ext != ".tck") - outFile += ".fib"; - mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp index 8e2b643adf..e8c3234e99 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.cpp @@ -1,264 +1,230 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkDenoisingView.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const std::string QmitkDenoisingView::VIEW_ID = "org.mitk.views.denoisingview"; QmitkDenoisingView::QmitkDenoisingView() : QmitkAbstractView() , m_Controls( 0 ) , m_ImageNode(nullptr) - , m_SelectedFilter(TV) { } QmitkDenoisingView::~QmitkDenoisingView() { } void QmitkDenoisingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkDenoisingViewControls; m_Controls->setupUi( parent ); - m_Controls->m_SelectFilterComboBox->clear(); - m_Controls->m_SelectFilterComboBox->insertItem(TV, "Total-variation filter"); - m_Controls->m_SelectFilterComboBox->insertItem(GAUSS, "Discrete gaussian filter"); - m_Controls->m_SelectFilterComboBox->insertItem(NLM, "Non-local means filter"); -// m_Controls->m_SelectFilterComboBox->insertItem(NLM_MORITZ, "Non-local means filter"); - connect( (QObject*)(m_Controls->m_ApplyButton), SIGNAL(clicked()), this, SLOT(StartDenoising())); - connect( (QObject*)(m_Controls->m_SelectFilterComboBox), SIGNAL(activated(int)), this, SLOT(UpdateGui(int))); - connect( (QObject*)(m_Controls->m_InputImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui(int))); + connect( (QObject*)(m_Controls->m_SelectFilterComboBox), SIGNAL(activated(int)), this, SLOT(UpdateGui())); + connect( (QObject*)(m_Controls->m_InputImageBox), SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui())); m_Controls->m_InputImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); m_Controls->m_InputImageBox->SetPredicate( isImagePredicate ); - UpdateGui(0); + UpdateGui(); } } void QmitkDenoisingView::SetFocus() { m_Controls->m_SelectFilterComboBox->setFocus(); } void QmitkDenoisingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& ) { } void QmitkDenoisingView::StartDenoising() { typedef itk::NonLocalMeansDenoisingFilter< DiffusionPixelType > NlmRicianFilterType; typedef itk::DiscreteGaussianImageFilter < DwiVolumeType, DwiVolumeType > GaussianFilterType; typedef itk::PatchBasedDenoisingImageFilter < DwiVolumeType, DwiVolumeType > NlmFilterType; typedef itk::VectorImageToImageFilter < DiffusionPixelType > ExtractFilterType; typedef itk::ComposeImageFilter < itk::Image > ComposeFilterType; m_ImageNode = m_Controls->m_InputImageBox->GetSelectedNode(); mitk::Image::Pointer input_image = dynamic_cast(m_ImageNode->GetData()); if (!mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(input_image)) return; DwiImageType::Pointer itkVectorImagePointer = mitk::DiffusionPropertyHelper::GetItkVectorImage(input_image); mitk::Image::Pointer denoised_image = nullptr; mitk::DataNode::Pointer denoised_image_node = mitk::DataNode::New(); - switch (m_SelectedFilter) + switch (m_Controls->m_SelectFilterComboBox->currentIndex()) { - case TV: + case 0: { ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i=0; iGetVectorLength(); ++i) { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); extractor->SetIndex( i ); extractor->Update(); DwiVolumeType::Pointer gradient_volume = extractor->GetOutput(); itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::Pointer filter = itk::TotalVariationDenoisingImageFilter< DwiVolumeType, DwiVolumeType >::New(); filter->SetInput(gradient_volume); filter->SetLambda(m_Controls->m_TvLambdaBox->value()); filter->SetNumberIterations(m_Controls->m_TvIterationsBox->value()); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); denoised_image_node->SetName( "TotalVariation" ); break; } - case GAUSS: + case 1: { ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); GaussianFilterType::Pointer filter = GaussianFilterType::New(); filter->SetVariance(m_Controls->m_GaussVarianceBox->value()); filter->SetInput(extractor->GetOutput()); filter->Update(); composer->SetInput(i, filter->GetOutput()); } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); denoised_image_node->SetName( "Gauss" ); break; } - case NLM: + case 2: { typedef itk::Statistics::GaussianRandomSpatialNeighborSubsampler< NlmFilterType::PatchSampleType, DwiVolumeType::RegionType > SamplerType; // sampling the image to find similar patches SamplerType::Pointer sampler = SamplerType::New(); sampler->SetRadius( m_Controls->m_NlmSearchRadiusBox->value() ); sampler->SetVariance( m_Controls->m_NlmSearchRadiusBox->value()*m_Controls->m_NlmSearchRadiusBox->value() ); sampler->SetNumberOfResultsRequested( m_Controls->m_NlmNumPatchesBox->value() ); MITK_INFO << "Starting NLM denoising"; ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput( itkVectorImagePointer ); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < itkVectorImagePointer->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); NlmFilterType::Pointer filter = NlmFilterType::New(); filter->SetInput(extractor->GetOutput()); filter->SetPatchRadius(m_Controls->m_NlmPatchSizeBox->value()); filter->SetNoiseModel(NlmFilterType::RICIAN); filter->UseSmoothDiscPatchWeightsOn(); filter->UseFastTensorComputationsOn(); filter->SetNumberOfIterations(m_Controls->m_NlmIterationsBox->value()); filter->SetSmoothingWeight( 1 ); filter->SetKernelBandwidthEstimation(true); filter->SetSampler( sampler ); filter->Update(); composer->SetInput(i, filter->GetOutput()); MITK_INFO << "Gradient " << i << " finished"; } composer->Update(); denoised_image = mitk::GrabItkImageMemory(composer->GetOutput()); denoised_image_node->SetName( "NLM" ); break; } - case NLM_MORITZ: - { - NlmRicianFilterType::Pointer nlmFilter = NlmRicianFilterType::New(); - nlmFilter->SetInputImage( itkVectorImagePointer ); -// nlmFilter->SetUseRicianAdaption(m_Controls->m_RicianCheckbox->isChecked()); -// nlmFilter->SetUseJointInformation(m_Controls->m_JointInformationCheckbox->isChecked()); - nlmFilter->SetSearchRadius(m_Controls->m_NlmSearchRadiusBox->value()); - nlmFilter->SetComparisonRadius(m_Controls->m_NlmPatchSizeBox->value()); -// nlmFilter->SetVariance(m_Controls->m_NlmVarianceBox->value()); - nlmFilter->Update(); - denoised_image = mitk::GrabItkImageMemory(nlmFilter->GetOutput()); - denoised_image_node->SetName( "NLM-Rician" ); - break; - } } denoised_image->GetPropertyList()->ReplaceProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( static_cast( input_image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ) ); denoised_image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(input_image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); mitk::DiffusionPropertyHelper propertyHelper( denoised_image ); propertyHelper.InitializeImage(); denoised_image_node->SetData( denoised_image ); GetDataStorage()->Add(denoised_image_node, m_ImageNode); } -void QmitkDenoisingView::UpdateGui(int filter) +void QmitkDenoisingView::UpdateGui() { m_Controls->m_ApplyButton->setEnabled(false); m_Controls->m_TvFrame->setVisible(false); m_Controls->m_NlmFrame->setVisible(false); m_Controls->m_GaussFrame->setVisible(false); - switch (filter) + switch (m_Controls->m_SelectFilterComboBox->currentIndex()) { case 0: { - m_SelectedFilter = TV; m_Controls->m_TvFrame->setVisible(true); break; } case 1: { - m_SelectedFilter = GAUSS; m_Controls->m_GaussFrame->setVisible(true); break; } case 2: { - m_SelectedFilter = NLM; m_Controls->m_NlmFrame->setVisible(true); break; } - case 3: - { - m_SelectedFilter = NLM_MORITZ; - m_Controls->m_NlmFrame->setVisible(true); - break; - } - default : - { - m_SelectedFilter = TV; - } } if (m_Controls->m_InputImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_ApplyButton->setEnabled(true); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.h index 366275705d..066c6f3aae 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingView.h @@ -1,90 +1,78 @@ /*=================================================================== 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 _QMITKQmitkDenoisingView_H_INCLUDED #define _QMITKQmitkDenoisingView_H_INCLUDED #include #include #include "ui_QmitkDenoisingViewControls.h" #include #include #include #include #include /** * \class QmitkDenoisingView * \brief View for denoising diffusion-weighted images. */ -class QmitkDenoisingView; - - class QmitkDenoisingView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkDenoisingView(); virtual ~QmitkDenoisingView(); /** Typedefs */ typedef short DiffusionPixelType; typedef itk::VectorImage DwiImageType; typedef itk::Image DwiVolumeType; virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; private slots: void StartDenoising(); - void UpdateGui(int filter); + void UpdateGui(); private: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkDenoisingViewControls* m_Controls; mitk::DataNode::Pointer m_ImageNode; - - enum FilterType { - TV, - GAUSS, - NLM, - NLM_MORITZ - } m_SelectedFilter; - - friend class QmitkDenoisingWorker; }; #endif // _QmitkDenoisingView_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingViewControls.ui index 35da1a192a..340b789ecd 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.denoising/src/internal/QmitkDenoisingViewControls.ui @@ -1,329 +1,345 @@ QmitkDenoisingViewControls 0 0 351 734 0 0 QmitkTemplate 6 9 9 9 9 QFrame::NoFrame QFrame::Raised 0 0 0 0 Input Image: - + + + + Total-variation + + + + + Discrete Gaussian + + + + + Non-local means + + + Filter: Parameters 0 0 0 0 Lambda: 4 0.000000000000000 1.000000000000000 0.100000000000000 0.100000000000000 Iterations: 1 999 1 Parameters 0 0 0 0 Sampling Radius: 1 999 10 1 999 4 1 999 4 Num. Patches: 1 50 1 Iterations: Patch Size: Parameters 0 0 0 0 4 1.000000000000000 99999.000000000000000 0.100000000000000 1.000000000000000 Variance: Start Qt::Vertical 20 40 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox index eb0ed3b82f..ca828cbd50 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox @@ -1,74 +1,76 @@ /** \page org_mitk_views_streamlinetracking Streamline Tractography This view enables streamline tractography on various input data. The corresponding command line application is named "MitkStreamlineTractography". Available sections: - \ref StrTrackUserManualInputData - \ref StrTrackUserManualInteractive - \ref StrTrackUserManualParameters - \ref StrTrackUserManualAdvancedParameters - \ref StrTrackUserManualReferences \section StrTrackUserManualInputData Input Data Input: Select the data you want to track on in the datamanager. Supported file types are: \li One or multiple DTI images selected in the datamanager. \li One ODF image, e.g. obtained using MITK Q-ball reconstruction or MRtrix CSD (tractography similar to [6]). \li One peak image (4D float image). \li One raw diffusion-weighted image for machine learning based tractography [1]. Optional Input: -\li Binary mask used to define the seed voxels. If no seed mask is specified, the whole image volume is seeded. -\li Binary mask used to constrain the generated streamlines. Streamlines can not leave the mask area. +\li Mask used to define the seed voxels. If no seed mask is specified, the whole image volume is seeded. +\li Mask used to constrain the generated streamlines. Streamlines can not leave the mask area. +\li Mask to define the target region. A streamline is only considered valid if it starts and ends in the target and seed region region. Without seed region, both endpoints need to be located in the target region. \li Binary mask used to define stopping regions. Streamlines that enter the mask area are stopped immediately. -\li Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. \li FA/GFA image used to determine streamline termination. If no image is specified, the FA/GFA image is automatically calculated from the input image. If multiple tensor images are used as input, it is recommended to provide such an image since the FA maps calculated from the individual input tensor images can not provide a suitable termination criterion. \li Tractography Forest: Needed for machine learning based tractography [1]. \section StrTrackUserManualInteractive Interactive Tractography Interactive tractography enables the dynamic placement of spherical seed regions simply by clicking into the image (similar to [5]). Parameters are the number of seed points and the radius of the spherical seed region. The seed points are randomly distributed inside the sphere around the selected position. By clicking and holding the left mouse button while moving the mouse around the image, the fiber connections originating from the selected region can be explored dynamically. When "Update on Parameter Change" is checked, each parameter change causes an instant retracking with the new parameters. This enables an intuitive exploration of the effcts that the individual parameters have on the resulting tractogram. \section StrTrackUserManualParameters Parameters \li Mode: Toggle between deterministic and probabilistic tractography. Peak tracking only supports deterministic mode. The probabilistic method simply samples the output direction from the discrete probability ditribution provided by the discretized ODF. \li Seeds per voxel: If set to 1, the seed is defined as the voxel center. If > 1 the seeds are distributet randomly inside the voxel. \li Max. num. fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. \li Cutoff: If the streamline reaches a position with an FA value or peak magnitude lower than the speciefied threshold, tracking is terminated. Typical values are 0.2 for FA/GFA and 0.1 for CSD peaks. \li ODF Cutoff: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. For MRtrix CSD fODF images, a typical value is 0.1. \li Sharpen ODFs: If you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). This is not necessary (and not recommended) for CSD fODFs, since they are naturally much sharper. \section StrTrackUserManualAdvancedParameters Advanced Parameters \li Step Size: The algorithm proceeds along the streamline with a fixed stepsize. Default is 0.5*minSpacing. \li Angular threshold: Maximum angle between two successive steps (in degree). Default is 90° * step_size. For probabilistic tractography, candidate directions exceeding this threshold have probability 0, i.e. the respective ODF value is set to zero. The probabilities of the valid directions are normalized to sum to 1. \li Min. Tract Length: Shorter fibers are discarded. +\li Loop Check: Stop streamline if the threshold on the angular stdev over the last 4 voxel lengths is exceeded. -1 = no loop check. +\li Trials Per Seed: Try each seed N times until a valid streamline is obtained (only for probabilistic tractography). \li f and g values to balance between FACT [2] and TEND [3,4] tracking (only for tensor based tractography). For further information please refer to [2,3] \li Flip directions: Internally flips progression directions. This might be necessary depending on the input data. \li Enable Trilinear Interpolation: By default the image values are interpolated. Keep in mind that in the noninterpolated case, the TEND term is only applied once per voxel. In the interpolated case the TEND term is applied at each integration step which results in much higher curvatures and has to be compensated by an according choice of f and g. \li Enable Gray Matter Seeding: Seeds are onyl placed inside of the gray matter. Needs tissue label image. \section StrTrackUserManualNeighbourhoodSampling Neighbourhood Sampling (for details see [1]) \li Neighborhood Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. \li Sampling Distance: Distance of the sampling positions from the current streamline position (in voxels). \li Use Only Frontal Samples: Only neighborhood samples in front of the current streamline position are considered. \li Use Stop-Votes: If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. \section StrTrackUserManualPostprocessing Output and Postprocessing \li Compress Fibers: Whole brain tractograms obtained with a small step size can contain billions of points. The tractograms can be compressed by removing points that do not really contribute to the fiber shape, such as many points on a straight line. An error threshold (in mm) can be defined to specify which points should be removed and which not. \li Output Probability Map: No streamline are generated. Instead, the tractography outputs a probability map that indicates the probability of a fiber to reach a voxel from the selected seed region. For this measure to be sensible, the number of seeds per voxel needs to be rather large. \section StrTrackUserManualReferences References [1] Neher, Peter F., Marc-Alexandre Côté, Jean-Christophe Houde, Maxime Descoteaux, and Klaus H. Maier-Hein. “Fiber Tractography Using Machine Learning.” NeuroImage. Accessed July 19, 2017. doi:10.1016/j.neuroimage.2017.07.028.\n [2] Mori, Susumu, Walter E. Kaufmann, Godfrey D. Pearlson, Barbara J. Crain, Bram Stieltjes, Meiyappan Solaiyappan, and Peter C. M. Van Zijl. “In Vivo Visualization of Human Neural Pathways by Magnetic Resonance Imaging.” Annals of Neurology 47 (2000): 412–414.\n [3] Weinstein, David, Gordon Kindlmann, and Eric Lundberg. “Tensorlines: Advection-Diffusion Based Propagation through Diffusion Tensor Fields.” In Proceedings of the Conference on Visualization’99: Celebrating Ten Years, 249–253, n.d.\n [4] Lazar, Mariana, David M. Weinstein, Jay S. Tsuruda, Khader M. Hasan, Konstantinos Arfanakis, M. Elizabeth Meyerand, Benham Badie, et al. “White Matter Tractography Using Diffusion Tensor Deflection.” Human Brain Mapping 18, no. 4 (2003): 306–321.\n [5] Chamberland, M., K. Whittingstall, D. Fortin, D. Mathieu, and M. Descoteaux. “Real-Time Multi-Peak Tractography for Instantaneous Connectivity Display.” Front Neuroinform 8 (2014): 59. doi:10.3389/fninf.2014.00059.\n [6] Tournier, J-Donald, Fernando Calamante, and Alan Connelly. “MRtrix: Diffusion Tractography in Crossing Fiber Regions.” International Journal of Imaging Systems and Technology 22, no. 1 (March 2012): 53–66. doi:10.1002/ima.22005. */ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index b33d3b55eb..57ef470c54 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,838 +1,850 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include #include // VTK #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) : m_View(view) { } void QmitkStreamlineTrackingWorker::run() { m_View->m_Tracker->Update(); m_View->m_TrackingThread.quit(); } QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_TrackingWorker(this) , m_Controls(nullptr) , m_FirstTensorProbRun(true) , m_FirstInteractiveRun(true) , m_TrackingHandler(nullptr) , m_ThreadIsRunning(false) , m_DeleteTrackingHandler(false) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TargetImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ForestBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestBox->SetPredicate(isTractographyForest); m_Controls->m_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageBox->SetZeroEntryText("--"); m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_SeedImageBox->SetZeroEntryText("--"); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_StopImageBox->SetZeroEntryText("--"); m_Controls->m_TargetImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_TargetImageBox->SetZeroEntryText("--"); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_TargetImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(valueChanged(double)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(valueChanged(int)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskInterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_LoopCheckBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); + connect( m_Controls->m_TrialsPerSeedBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); StartStopTrackingGui(false); } UpdateGui(); } void QmitkStreamlineTrackingView::StopTractography() { if (m_Tracker.IsNull()) return; m_Tracker->SetStopTracking(true); } void QmitkStreamlineTrackingView::TimerUpdate() { if (m_Tracker.IsNull()) return; QString status_text(m_Tracker->GetStatusText().c_str()); m_Controls->m_StatusTextBox->setText(status_text); } void QmitkStreamlineTrackingView::BeforeThread() { m_TrackingTimer->start(1000); } void QmitkStreamlineTrackingView::AfterThread() { m_TrackingTimer->stop(); if (!m_Tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_ParentNode->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked() && fiberBundle->GetNumberOfLines()>0) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); m_Tracker->SetDicomProperties(fib); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_ParentNode->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); GetDataStorage()->Add(node, m_ParentNode); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_ParentNode->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_ParentNode); } } if (m_InteractivePointSetNode.IsNotNull()) m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); StartStopTrackingGui(false); if (m_DeleteTrackingHandler) DeleteTrackingHandler(); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if (m_ThreadIsRunning) return; if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked())) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,0,0)); DoFiberTracking(); } void QmitkStreamlineTrackingView::OnParameterChanged() { if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { if (m_FirstInteractiveRun) { QMessageBox::information(nullptr, "Information", "Place and move a spherical seed region anywhere in the image by left-clicking and dragging. If the seed region is colored red, tracking is in progress. If the seed region is colored white, tracking is finished.\nPlacing the seed region for the first time in a newly selected dataset might cause a short delay, since the tracker needs to be initialized."); m_FirstInteractiveRun = false; } QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); connect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; m_SliceChangeListener.RenderWindowPartActivated(this->GetRenderWindowPart()); disconnect(&m_SliceChangeListener, SIGNAL(SliceChanged()), this, SLOT(OnSliceChanged())); } } void QmitkStreamlineTrackingView::OnSliceChanged() { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (!m_ThreadIsRunning && m_TrackingHandler != nullptr) { delete m_TrackingHandler; m_TrackingHandler = nullptr; m_DeleteTrackingHandler = false; } else if (m_ThreadIsRunning) { m_DeleteTrackingHandler = true; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { std::vector< mitk::DataNode::Pointer > last_nodes = m_InputImageNodes; m_InputImageNodes.clear(); m_InputImages.clear(); m_AdditionalInputImages.clear(); bool retrack = false; for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData())) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr) { int dim = img->GetDimension(); unsigned int* dimensions = img->GetDimensions(); if (dim==4 && dimensions[3]%3==0) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); retrack = true; } else if (dim==3) { m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } } } // sometimes the OnSelectionChanged event is sent twice and actually no selection has changed for the first event. We need to catch that. if (last_nodes.size() == m_InputImageNodes.size()) { bool same_nodes = true; for (unsigned int i=0; im_TensorImageLabel->setText("mandatory"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageBox->setEnabled(false); m_Controls->mFaImageLabel->setEnabled(false); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestBox->setEnabled(false); m_Controls->m_ForestLabel->setEnabled(false); m_Controls->commandLinkButton->setEnabled(false); + m_Controls->m_TrialsPerSeedBox->setEnabled(false); + m_Controls->m_TrialsPerSeedLabel->setEnabled(false); + + // trials per seed are only important for probabilistic tractography + if (m_Controls->m_ModeBox->currentIndex()==1) + { + m_Controls->m_TrialsPerSeedBox->setEnabled(true); + m_Controls->m_TrialsPerSeedLabel->setEnabled(true); + } if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText( ( std::to_string(m_InputImageNodes.size()) + " images selected").c_str() ); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->m_InputData->setTitle("Input Data"); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestBox->setEnabled(true); m_Controls->m_ForestLabel->setEnabled(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } else m_Controls->m_InputData->setTitle("Please Select Input Data"); } void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) { m_ThreadIsRunning = start; if (!m_Controls->m_InteractiveBox->isChecked()) { m_Controls->commandLinkButton_2->setVisible(start); m_Controls->commandLinkButton->setVisible(!start); m_Controls->m_InteractiveBox->setEnabled(!start); m_Controls->m_StatusTextBox->setVisible(start); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_ThreadIsRunning) return; if (m_InputImages.empty()) return; if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) return; StartStopTrackingGui(true); m_Tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { typedef mitk::ImageToItk CasterType; if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImages.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); StartStopTrackingGui(false); return; } if (m_FirstTensorProbRun) { QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. ODFs are sharpened (min-max normalized and raised to the power of 4). TEND parameters are ignored."); m_FirstTensorProbRun = false; } if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TensorImage::ItkTensorImageType::Pointer itkImg = mitk::TensorImage::ItkTensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); typedef itk::TensorImageToOdfImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkImg ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); dynamic_cast(m_TrackingHandler)->SetIsOdfFromTensor(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (int i=0; i<(int)m_InputImages.size(); i++) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(i)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::ConstPointer itkImg = caster->GetOutput(); dynamic_cast(m_TrackingHandler)->AddTensorImage(itkImg); } if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(m_TrackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); dynamic_cast(m_TrackingHandler)->SetOdfImage(itkImg); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestBox->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography (raw dMRI tractography) selected. Did you accidentally select the raw diffusion-weighted image in the datamanager?"); StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestBox->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { int num_previous_directions = (forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { int num_previous_directions = (forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); StartStopTrackingGui(false); return; } } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); StartStopTrackingGui(false); return; } try { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(0)); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(itkImg); } dynamic_cast(m_TrackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); } catch(...) { QMessageBox::information(nullptr, "Error", "Peak tracker could not be initialized. Is your input image in the correct format (4D float image, peaks in the 4th dimension)?"); StartStopTrackingGui(false); return; } } m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { m_Tracker->SetSeedPoints(m_SeedPoints); } else if (m_Controls->m_SeedImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetSeedImage(mask); } m_Tracker->SetInterpolateMask(false); if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetMaskImage(mask); m_Tracker->SetInterpolateMask(m_Controls->m_MaskInterpolationBox->isChecked()); } if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetStoppingRegions(mask); } if (m_Controls->m_TargetImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetTargetRegions(mask); } m_Tracker->SetVerbose(!m_Controls->m_InteractiveBox->isChecked()); m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); - m_Tracker->SetAposterioriCurvCheck(false); + m_Tracker->SetTrialsPerSeed(m_Controls->m_TrialsPerSeedBox->value()); m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); + m_Tracker->SetLoopCheck(m_Controls->m_LoopCheckBox->value()); m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); m_ParentNode = m_InputImageNodes.at(0); m_TrackingThread.start(QThread::LowestPriority); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui index ff66c862a0..f5a161760e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui @@ -1,1131 +1,1189 @@ QmitkStreamlineTrackingViewControls 0 0 413 1385 0 0 QmitkTemplate 0 0 3 3 Please Select Input Data 0 0 0 0 Random forest for machine learning based tractography. QComboBox::AdjustToMinimumContentsLength - Stop Image: Fibers that enter a region defined in this image will stop immediately. QComboBox::AdjustToMinimumContentsLength - Input Image. ODF, tensor, peak, and, in case of ML tractography, raw diffusion-weighted images are currently supported. <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. QComboBox::AdjustToMinimumContentsLength - Mask Image: Input Image. ODF, tensor and peak images are currently supported. Input Image: If an image is selected, the stopping criterion is not calculated from the input image but instead the selected image is used. QComboBox::AdjustToMinimumContentsLength - Target Image: FA/GFA Image: - Fibers that do not start and end in a region defined in this image will be discarded. + Fibers that do not start and end in a region defined in this image or the seed image will be discarded. QComboBox::AdjustToMinimumContentsLength - Seed Image: Seed points are only placed inside the mask image. If no seed mask is selected, the whole image is seeded. QComboBox::AdjustToMinimumContentsLength - Tractography Forest: false Start Tractography 0 0 Interactive Tractography 0 0 0 0 Number of seed points normally distributed around selected position. 1 9999999 50 Num. Seeds: true Dynamically pick seed location by click into image. Enable Interactive Tractography Seedpoints are normally distributed within a sphere centered at the selected position with the specified radius (in mm). 2 50.000000000000000 0.100000000000000 2.500000000000000 Radius: true When checked, parameter changes cause instant retracking while in interactive mode. Update on Parameter Change true 0 0 Parameters 0 0 0 0 Mode: Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. Deterministic Probabilistic Seeds per Voxel: Number of seed points placed in each voxel. 1 9999999 Max. num. fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. -1 999999999 -1 Cutoff: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 ODF Cutoff: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. 5 1.000000000000000 0.100000000000000 0.100000000000000 If you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). This is not necessary for CSD fODFs, since they are naturally much sharper. Sharpen ODFs Qt::Horizontal QSizePolicy::Fixed 200 0 Advanced Parameters 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 - - - - - - - Min. Tract Length: - - - - - + + - Step size (in voxels) + f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 - - 0.010000000000000 - - 10.000000000000000 + 1.000000000000000 0.100000000000000 - 0.500000000000000 + 1.000000000000000 - - + + - + Default: 90° * step_size - - Angular Threshold: + + -1 - - - - - - + + 90 - - Flip directions: + + 1 + + + -1 - - + + - f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). + Try each seed N times until a valid streamline is obtained (only for probabilistic). - - 2 + + Minimum fiber length (in mm) - - 1.000000000000000 + + 1 - - 0.100000000000000 + + 999 - 1.000000000000000 + 10 - + f parameter of tensor tractography. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). f: - - - - - - - g: - - - - - - - Shorter fibers are discarded. - - - Minimum fiber length (in mm) - - - 1 - - - 999.000000000000000 - - - 1.000000000000000 - - - 20.000000000000000 - - - - + QFrame::NoFrame QFrame::Raised 0 0 0 0 Internally flips progression directions. This might be necessary depending on the input data. x Internally flips progression directions. This might be necessary depending on the input data. y Internally flips progression directions. This might be necessary depending on the input data. z - - + + - Step Size: + Flip directions: - - + + - Default: 90° * step_size + - - -1 + + Angular Threshold: + + + + + + + f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). + + + 2 - 90 + 1.000000000000000 - 1 + 0.100000000000000 - -1 + 0.000000000000000 - - + + - f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). + Step size (in voxels) 2 + + 0.010000000000000 + - 1.000000000000000 + 10.000000000000000 0.100000000000000 - 0.000000000000000 + 0.500000000000000 + + + + + + + + + + Min. Tract Length: + + + + + + + + + + Step Size: + + + + + + + + + + Trials Per Seed: + + + + + + + Shorter fibers are discarded. + + + Minimum fiber length (in mm) + + + 1 + + + 999.000000000000000 + + + 1.000000000000000 + + + 20.000000000000000 + + + + + + + + + + g: + + + + + + + + + + Loop Check: + + + + + + + Shorter fibers are discarded. + + + Maximum allowed angular SDTEV over 4 voxel lengths. Default: no loop check. + + + -1 + + + 180 + + + -1 QFrame::NoFrame QFrame::Raised 0 0 0 0 If false, nearest neighbor interpolation is used. Interpolate Tractography Data true The tractography mask is not interpolated by default. Interpolate ROI Images true 0 0 Neighborhood Sampling 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Num. Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 Sampling Distance: Sampling distance (in voxels) 2 10.000000000000000 0.100000000000000 0.250000000000000 Only neighborhood samples in front of the current streamline position are considered. Use Only Frontal Samples false If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. Use Stop-Votes false Qt::Vertical QSizePolicy::Expanding 20 220 0 0 Output and Postprocessing 0 0 0 0 QFrame::NoFrame QFrame::Raised 0 0 0 0 Resample fibers using the specified error constraint. Compress Fibers true Qt::StrongFocus Lossy fiber compression. Recommended for large tractograms. Maximum error in mm. 3 10.000000000000000 0.010000000000000 0.100000000000000 Output probability map instead of tractogram. Output Probability Map false true Stop Tractography true 0 0 true QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h
 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h