diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 8fbef82e1f..1e5c3dd017 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1093 +1,906 @@ /*=================================================================== 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 #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_TissueImage(nullptr) , m_OutputProbabilityMap(nullptr) , 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_WmLabel(3) // mrtrix 5ttseg labels - , m_GmLabel(1) // mrtrix 5ttseg labels - , m_SeedOnlyGm(false) - , m_ControlGmEndings(false) , m_MaxNumTracts(-1) , m_Verbose(true) , m_AposterioriCurvCheck(false) , m_DemoMode(false) , m_Random(true) , m_UseOutputProbabilityMap(false) , m_CurrentTracts(0) , m_Progress(0) , m_StopTracking(false) - , m_InterpolateMask(false) + , m_InterpolateMask(true) { this->SetNumberOfRequiredInputs(0); - - m_MaskInterpolator = itk::LinearInterpolateImageFunction< itk::Image< unsigned char, 3 >, float >::New(); - m_StopInterpolator = itk::LinearInterpolateImageFunction< itk::Image< unsigned char, 3 >, float >::New(); } 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 = ItkUcharImgType::New(); + 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_TargetRegions = ItkUintImgType::New(); + 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_SeedImage = ItkUcharImgType::New(); + 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 = ItkUcharImgType::New(); + 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(); - else - m_SeedOnlyGm = false; - - if (m_TissueImage.IsNull()) - { - if (m_SeedOnlyGm) - { - MITK_WARN << "StreamlineTracking - Cannot seed in gray matter. No tissue type image set."; - m_SeedOnlyGm = false; - } - - if (m_ControlGmEndings) - { - MITK_WARN << "StreamlineTracking - Cannot control gray matter endings. No tissue type image set."; - m_ControlGmEndings = false; - } - } - else - { - if (m_ControlGmEndings) - m_SeedOnlyGm = true; - if (m_ControlGmEndings || m_SeedOnlyGm) - std::cout << "StreamlineTracking - Using tissue image" << std::endl; - else - MITK_WARN << "StreamlineTracking - Tissue image set but no gray matter seeding or fiber endpoint-control enabled" << std::endl; - } 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_SeedOnlyGm && m_ControlGmEndings) - InitGrayMatterEndings(); - 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; 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 - 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 - 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 - 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::InitGrayMatterEndings() -{ - m_TrackingHandler->SetAngularThreshold(0); - m_GmStubs.clear(); - if (m_TissueImage.IsNotNull()) - { - std::cout << "StreamlineTracking - initializing GM endings" << std::endl; - ImageRegionConstIterator< ItkUcharImgType > it(m_TissueImage, m_TissueImage->GetLargestPossibleRegion() ); - it.GoToBegin(); - - vnl_vector_fixed d1; d1.fill(0.0); - std::deque< vnl_vector_fixed > olddirs; - while (olddirs.size() start; - m_TissueImage->TransformIndexToPhysicalPoint(s_idx, start); - itk::Point wm_p; - float max = -1; - FiberType fib; - - for (int x : {-1,0,1}) - for (int y : {-1,0,1}) - for (int z : {-1,0,1}) - { - if (x==y && y==z) - continue; - - ItkUcharImgType::IndexType e_idx; - e_idx[0] = s_idx[0] + x; - e_idx[1] = s_idx[1] + y; - e_idx[2] = s_idx[2] + z; - - if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_WmLabel ) - continue; - - itk::ContinuousIndex end; - m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); - - d1 = m_TrackingHandler->ProposeDirection(end, olddirs, s_idx); - if (d1.magnitude()<0.0001) - continue; - d1.normalize(); - - vnl_vector_fixed< float, 3 > d2; - d2[0] = end[0] - start[0]; - d2[1] = end[1] - start[1]; - d2[2] = end[2] - start[2]; - d2.normalize(); - float a = fabs(dot_product(d1,d2)); - if (a>max) - { - max = a; - wm_p = end; - } - } - - if (max>=0) - { - fib.push_back(start); - fib.push_back(wm_p); - m_GmStubs.push_back(fib); - } - } - ++it; - } - } - m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); -} - - 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; } -bool StreamlineTrackingFilter -::IsInGm(const itk::Point &pos) -{ - if (m_TissueImage.IsNull()) - return true; - - ItkUcharImgType::IndexType idx; - m_TissueImage->TransformPhysicalPointToIndex(pos, idx); - if (m_TissueImage->GetLargestPossibleRegion().IsInside(idx) && m_TissueImage->GetPixel(idx) == m_GmLabel) - return true; - - return false; -} - 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)) + 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)) + 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)) + 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) { 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; } else // if yes, add new point to streamline { tractLength += m_StepSize; if (front) fib->push_front(pos); else fib->push_back(pos); 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 (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 m_Distance = 5; if (fib->size()<3) return 0; 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; 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(); vectors.push_back(v); if (c==0) meanV += v; c++; } } else { int c=fib->size()-1; while(dist0) { 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(); vectors.push_back(v); if (c==(int)fib->size()-1) 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(); } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - Calculating seed points."; m_SeedPoints.clear(); - if (!m_ControlGmEndings) - { - typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; - MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); - sit.GoToBegin(); + typedef ImageRegionConstIterator< ItkFloatImgType > MaskIteratorType; + MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); + sit.GoToBegin(); - while (!sit.IsAtEnd()) + while (!sit.IsAtEnd()) + { + if (sit.Value()>0) { - 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) ) { - ItkUcharImgType::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_SeedOnlyGm || IsInGm(worldPos))) + m_SeedPoints.push_back(worldPos); + for (int s = 1; s < m_SeedsPerVoxel; s++) { - 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); + 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); - } + itk::Point worldPos; + m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); + m_SeedPoints.push_back(worldPos); } } - ++sit; } - } - else - { - for (auto s : m_GmStubs) - m_SeedPoints.push_back(s[1]); + ++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)) + /// START DIR + // vnl_vector_fixed< float, 3 > 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) { - if (m_ControlGmEndings) - { - float a = dot_product(gm_start_dir, dir); - if (a<0) - dir = -dir; - } + /// 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); - if (m_ControlGmEndings) - { - fib.push_front(m_GmStubs[temp_i][0]); - CheckFiberForGmEnding(&fib); - } - else - { - // backward tracking (only if we don't explicitely start in the GM) - tractLength = FollowStreamline(worldPos, -dir, &fib, tractLength, true); - if (m_ControlGmEndings) - { - CheckFiberForGmEnding(&fib); - std::reverse(fib.begin(),fib.end()); - CheckFiberForGmEnding(&fib); - } - } + // backward tracking (only if we don't explicitely start in the GM) + tractLength = FollowStreamline(worldPos, -dir, &fib, tractLength, true); + counter = fib.size(); if (tractLength>=m_MinTractLength && counter>=2) { - ItkUintImgType::IndexType idx_begin, idx_end; - m_TargetRegions->TransformPhysicalPointToIndex(fib.front(), idx_begin); - m_TargetRegions->TransformPhysicalPointToIndex(fib.back(), idx_end); #pragma omp critical - if ( !m_TargetRegions->GetLargestPossibleRegion().IsInside(idx_end) || - !m_TargetRegions->GetLargestPossibleRegion().IsInside(idx_begin) || - (m_TargetRegions->GetPixel(idx_begin)>0 && m_TargetRegions->GetPixel(idx_end)==m_TargetRegions->GetPixel(idx_begin)) ) + if ( IsValidFiber(&fib) ) { 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)) { if (!m_StopTracking) { std::cout << " \r"; MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; } m_StopTracking = true; } } } } } this->AfterTracking(); } -void StreamlineTrackingFilter::CheckFiberForGmEnding(FiberType* fib) +bool StreamlineTrackingFilter::IsValidFiber(FiberType* fib) { - if (m_TissueImage.IsNull()) - return; - - // first check if the current fibe rendpoint is located inside of the white matter - // if not, remove last fiber point and repeat - bool in_wm = false; - while (!in_wm && fib->size()>2) + if (m_TargetImageSet && m_SeedImageSet) { - ItkUcharImgType::IndexType idx; - m_TissueImage->TransformPhysicalPointToIndex(fib->back(), idx); - if (m_TissueImage->GetPixel(idx)==m_WmLabel) - in_wm = true; - else - fib->pop_back(); + 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; } - if (fib->size()<3 || !in_wm) + else if (m_TargetImageSet) { - fib->clear(); - return; + if ( mitk::imv::IsInsideMask(fib->front(), m_InterpolateMask, m_TargetInterpolator) + && mitk::imv::IsInsideMask(fib->back(), m_InterpolateMask, m_TargetInterpolator) ) + return true; + return false; } - - // get fiber direction at end point - vnl_vector_fixed< float, 3 > d1; - d1[0] = fib->back()[0] - fib->at(fib->size()-2)[0]; - d1[1] = fib->back()[1] - fib->at(fib->size()-2)[1]; - d1[2] = fib->back()[2] - fib->at(fib->size()-2)[2]; - d1.normalize(); - - // find closest gray matter voxel - ItkUcharImgType::IndexType s_idx; - m_TissueImage->TransformPhysicalPointToIndex(fib->back(), s_idx); - itk::Point gm_endp; - float max = -1; - - for (int x : {-1,0,1}) - for (int y : {-1,0,1}) - for (int z : {-1,0,1}) - { - if (x==y && y==z) - continue; - - ItkUcharImgType::IndexType e_idx; - e_idx[0] = s_idx[0] + x; - e_idx[1] = s_idx[1] + y; - e_idx[2] = s_idx[2] + z; - - if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_GmLabel ) - continue; - - itk::ContinuousIndex end; - m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); - vnl_vector_fixed< float, 3 > d2; - d2[0] = end[0] - fib->back()[0]; - d2[1] = end[1] - fib->back()[1]; - d2[2] = end[2] - fib->back()[2]; - d2.normalize(); - float a = dot_product(d1,d2); - if (a>max) - { - max = a; - gm_endp = end; - } - } - - if (max>=0) - fib->push_back(gm_endp); - else // no gray matter enpoint found -> delete fiber - fib->clear(); + 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 68611d17c7..8fc5b836a2 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,231 +1,224 @@ /*=================================================================== 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< 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, ItkUcharImgType::Pointer) ///< Seeds are only placed inside of this mask. - itkSetMacro( MaskImage, ItkUcharImgType::Pointer) ///< Tracking is only performed inside of this mask image. - itkSetMacro( TissueImage, ItkUcharImgType::Pointer) ///< + 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, ItkUcharImgType::Pointer) ///< Streamlines entering a stopping region will stop immediately - itkSetMacro( TargetRegions, ItkUintImgType::Pointer) ///< Only streamline starting and ending in this mask are retained + 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( SeedOnlyGm, bool ) ///< place seed points only in the gray matter - itkSetMacro( ControlGmEndings, 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( 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 ) ///< 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() {} - void InitGrayMatterEndings(); - void CheckFiberForGmEnding(FiberType* fib); + 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. - bool IsInGm(const itk::Point &pos); 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; - ItkUcharImgType::Pointer m_StoppingRegions; - ItkUintImgType::Pointer m_TargetRegions; - ItkUcharImgType::Pointer m_SeedImage; - ItkUcharImgType::Pointer m_MaskImage; - ItkUcharImgType::Pointer m_TissueImage; + ItkFloatImgType::Pointer m_StoppingRegions; + ItkFloatImgType::Pointer m_TargetRegions; + ItkFloatImgType::Pointer m_SeedImage; + ItkFloatImgType::Pointer m_MaskImage; ItkDoubleImgType::Pointer m_OutputProbabilityMap; 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_WmLabel; - int m_GmLabel; - bool m_SeedOnlyGm; - bool m_ControlGmEndings; int m_MaxNumTracts; bool m_Verbose; bool m_AposterioriCurvCheck; 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; void BuildFibers(bool check); int CheckCurvature(FiberType* 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< itk::Image< unsigned char, 3 >, float >::Pointer m_MaskInterpolator; - itk::LinearInterpolateImageFunction< itk::Image< unsigned char, 3 >, float >::Pointer m_StopInterpolator; + 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/Testing/mitkMachineLearningTrackingTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp index 25dd661f25..36f231a281 100644 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkMachineLearningTrackingTest.cpp @@ -1,108 +1,109 @@ /*=================================================================== 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 ItkUcharImgType; + 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; - ItkUcharImgType::Pointer seed; + 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 = ItkUcharImgType::New(); + 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 789e6d95d7..9e23f45457 100755 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTractographyTest.cpp @@ -1,430 +1,427 @@ /*=================================================================== 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 ItkUcharImgType; 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; - ItkUcharImgType::Pointer itk_seed_image; - ItkUcharImgType::Pointer itk_mask_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 = ItkUcharImgType::New(); + itk_seed_image = ItkFloatImgType::New(); mitk::CastToItkImage(seed_image, itk_seed_image); - itk_mask_image = ItkUcharImgType::New(); + 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->SetTissueImage(nullptr); - tracker->SetSeedOnlyGm(false); - tracker->SetControlGmEndings(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/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp index baa8cdfa0d..ce317deddf 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp @@ -1,485 +1,456 @@ /*=================================================================== 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("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("tissue_image", "", mitkCommandLineParser::String, "Tissue type image:", "image with tissue type labels (WM=3, GM=1)", 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("seed_gm", "", mitkCommandLineParser::Bool, "Seed only GM:", "Seed only in gray matter (requires tissue type image --tissue_image)"); - parser.addArgument("control_gm_endings", "", mitkCommandLineParser::Bool, "Control GM endings:", "Seed perpendicular to gray matter and enforce endings inside of the gray matter (requires tissue type image --tissue_image)"); parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded.", -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 seed_gm = false; - if (parsedArgs.count("seed_gm")) - seed_gm = us::any_cast(parsedArgs["seed_gm"]); - - bool control_gm_endings = false; - if (parsedArgs.count("control_gm_endings")) - control_gm_endings = us::any_cast(parsedArgs["control_gm_endings"]); - 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"]); 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"]); 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"]); - std::string tissueFile = ""; - if (parsedArgs.count("tissue_image")) - tissueFile = us::any_cast(parsedArgs["tissue_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"]); 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 mitkCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); - typedef itk::Image ItkUcharImgType; - typedef itk::Image ItkUintImgType; + 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); } - ItkUcharImgType::Pointer mask; + ItkFloatImgType::Pointer mask; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(maskFile)[0].GetPointer()); - mask = ItkUcharImgType::New(); + mask = ItkFloatImgType::New(); mitk::CastToItkImage(img, mask); } - ItkUcharImgType::Pointer seed; + ItkFloatImgType::Pointer seed; if (!seedFile.empty()) { MITK_INFO << "loading seed image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(seedFile)[0].GetPointer()); - seed = ItkUcharImgType::New(); + seed = ItkFloatImgType::New(); mitk::CastToItkImage(img, seed); } - ItkUcharImgType::Pointer stop; + ItkFloatImgType::Pointer stop; if (!stopFile.empty()) { MITK_INFO << "loading stop image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(stopFile)[0].GetPointer()); - stop = ItkUcharImgType::New(); + stop = ItkFloatImgType::New(); mitk::CastToItkImage(img, stop); } - ItkUintImgType::Pointer target; + ItkFloatImgType::Pointer target; if (!targetFile.empty()) { MITK_INFO << "loading target image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(targetFile)[0].GetPointer()); - target = ItkUintImgType::New(); + target = ItkFloatImgType::New(); mitk::CastToItkImage(img, target); } - ItkUcharImgType::Pointer tissue; - if (!tissueFile.empty()) - { - MITK_INFO << "loading tissue image"; - mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::Load(tissueFile)[0].GetPointer()); - tissue = ItkUcharImgType::New(); - mitk::CastToItkImage(img, tissue); - } - MITK_INFO << "loading additional images"; - typedef itk::Image ItkFloatImgType; 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->SetTissueImage(tissue); - tracker->SetSeedOnlyGm(seed_gm); - tracker->SetControlGmEndings(control_gm_endings); tracker->SetMaxNumTracts(max_tracts); tracker->SetTrackingHandler(handler); tracker->SetUseOutputProbabilityMap(output_prob_map); tracker->SetMinTractLength(min_tract_length); tracker->Update(); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (!output_prob_map) { 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.tractography/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp index be53ab681c..b33d3b55eb 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,859 +1,838 @@ /*=================================================================== 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_TissueImageBox->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(isBinaryPredicate, dimensionPredicate) ); + m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_SeedImageBox->SetZeroEntryText("--"); - m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); + m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); - m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); + m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_StopImageBox->SetZeroEntryText("--"); - m_Controls->m_TargetImageBox->SetPredicate( dimensionPredicate ); + m_Controls->m_TargetImageBox->SetPredicate( mitk::NodePredicateAnd::New(isImagePredicate, dimensionPredicate) ); m_Controls->m_TargetImageBox->SetZeroEntryText("--"); - m_Controls->m_TissueImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); - m_Controls->m_TissueImageBox->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_TissueImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); 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_SeedGmBox, 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()) ); 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_SeedGmBox->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()); - m_Controls->m_TissueImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); - m_Controls->label_10->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); - if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) - m_Controls->m_SeedGmBox->setEnabled(true); - else - m_Controls->m_SeedGmBox->setEnabled(false); - 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()) { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + 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()) { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + 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()) { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); + 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()) { - ItkUintImgType::Pointer mask = ItkUintImgType::New(); + ItkFloatImageType::Pointer mask = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TargetImageBox->GetSelectedNode()->GetData()), mask); m_Tracker->SetTargetRegions(mask); } - if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) - { - ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); - mitk::CastToItkImage(dynamic_cast(m_Controls->m_TissueImageBox->GetSelectedNode()->GetData()), mask); - m_Tracker->SetTissueImage(mask); - m_Tracker->SetSeedOnlyGm(m_Controls->m_SeedGmBox->isChecked()); - } - 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->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); m_Tracker->SetTrackingHandler(m_TrackingHandler); 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 83ee13c09d..ff66c862a0 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,1169 +1,1131 @@ QmitkStreamlineTrackingViewControls 0 0 413 1385 0 0 QmitkTemplate 0 0 3 3 Please Select Input Data 0 0 0 0 - - + + - Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. + Random forest for machine learning based tractography. QComboBox::AdjustToMinimumContentsLength - - - + + - Tissue Image: + Stop Image: - - + + - + Fibers that enter a region defined in this image will stop immediately. - - Tractography Forest: + + QComboBox::AdjustToMinimumContentsLength + + + - + + - - + + - + Input Image. ODF, tensor, peak, and, in case of ML tractography, raw diffusion-weighted images are currently supported. - Seed Image: - - - - - - - Input Image. ODF, tensor and peak images are currently supported. + <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> - - Input Image: + + true Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. QComboBox::AdjustToMinimumContentsLength - - - - - Fibers that enter a region defined in this image will stop immediately. - - - QComboBox::AdjustToMinimumContentsLength - - - - - - - - - - - + + - FA/GFA Image: + Mask Image: - - + + - Random forest for machine learning based tractography. + Input Image. ODF, tensor and peak images are currently supported. - - QComboBox::AdjustToMinimumContentsLength + + 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 - - - + + - Stop Image: + Target Image: - - + + - Mask Image: + FA/GFA Image: - - + + - Seed points are only placed inside the mask image. If no seed mask is selected, the whole image is seeded. + Fibers that do not start and end in a region defined in this image will be discarded. 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 - - - - - + + - Target Image: + Seed Image: - - + + - Fibers that do not start and end in a region defined in this image will be discarded. + 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) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Angular Threshold: Flip directions: f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 1.000000000000000 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: Default: 90° * step_size -1 90 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!). 2 1.000000000000000 0.100000000000000 0.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 - - - - Requires tissue image. - - - Enable Gray Matter Seeding - - - false - - - If false, nearest neighbor interpolation is used. Interpolate Tractography Data true The tractography mask is not interpolated by default. - Interpolate Mask Image + Interpolate ROI Images - false + 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