diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 6d578d7472..2c1475444e 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,943 +1,979 @@ /*=================================================================== 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" #define _USE_MATH_DEFINES #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_AngularThreshold(-1) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) , m_RandomSampling(false) , m_SamplingDistance(-1) , m_DeflectionMod(1.0) , m_OnlyForwardSamples(true) , m_UseStopVotes(true) , m_NumberOfSamples(30) , m_NumPreviousDirections(1) , m_StoppingRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_AposterioriCurvCheck(false) , m_AvoidStop(true) , m_DemoMode(false) , m_SeedOnlyGm(false) , m_ControlGmEndings(false) , m_WmLabel(3) // mrtrix 5ttseg labels , m_GmLabel(1) // mrtrix 5ttseg labels , m_StepSizeVox(-1) , m_SamplingDistanceVox(-1) , m_AngularThresholdDeg(-1) , m_MaxNumTracts(-1) , m_Random(true) , m_Verbose(true) + , m_UseOutputProbabilityMap(true) { this->SetNumberOfRequiredInputs(0); } void StreamlineTrackingFilter::BeforeTracking() { std::cout.setf(std::ios::boolalpha); 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); + } + if (m_StoppingRegions.IsNull()) { m_StoppingRegions = ItkUcharImgType::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; if (m_SeedImage.IsNull()) { m_SeedImage = ItkUcharImgType::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 std::cout << "StreamlineTracking - Using seed image" << std::endl; if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkUcharImgType::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; 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 ::IsValidPosition(itk::Point &pos) { ItkUcharImgType::IndexType idx; m_MaskImage->TransformPhysicalPointToIndex(pos, idx); if (!m_MaskImage->GetLargestPossibleRegion().IsInside(idx) || m_MaskImage->GetPixel(idx)==0) return false; return true; } float StreamlineTrackingFilter::GetRandDouble(float min, float max) { return (float)(rand()%((int)(10000*(max-min))) + 10000*min)/10000; } 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); ItkUcharImgType::IndexType idx; m_MaskImage->TransformPhysicalPointToIndex(pos, idx); if (IsValidPosition(pos)) { if (m_StoppingRegions->GetPixel(idx)>0) return direction; 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 (int i=0; i d; bool is_stop_voter = false; if (m_RandomSampling) { d[0] = GetRandDouble(); d[1] = GetRandDouble(); d[2] = GetRandDouble(); d.normalize(); d *= 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 (IsValidPosition(sample_pos)) 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 (IsValidPosition(sample_pos)) 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; iTransformPhysicalPointToIndex(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; } #pragma omp critical - if (m_DemoMode) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? + if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? { 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==fib->size()-1) meanV += v; c--; } } meanV.normalize(); for (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::GenerateData() { this->BeforeTracking(); std::vector< itk::Point > seedpoints; if (!m_ControlGmEndings) { typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion() ); MaskIteratorType mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); sit.GoToBegin(); mit.GoToBegin(); while( !sit.IsAtEnd() ) { if (sit.Value()==0 || mit.Value()==0 || (m_SeedOnlyGm && m_TissueImage->GetPixel(sit.GetIndex())!=m_GmLabel)) { ++sit; ++mit; continue; } for (int s=0; s start; if (m_SeedsPerVoxel>1) { start[0] = index[0]+GetRandDouble(-0.5, 0.5); start[1] = index[1]+GetRandDouble(-0.5, 0.5); start[2] = index[2]+GetRandDouble(-0.5, 0.5); } else { start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; } itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint( start, worldPos ); seedpoints.push_back(worldPos); } ++sit; ++mit; } } else { for (auto s : m_GmStubs) seedpoints.push_back(s[1]); } if (m_Random) { std::srand(std::time(0)); std::random_shuffle ( seedpoints.begin(), seedpoints.end() ); } bool stop = false; unsigned int current_tracts = 0; int num_seeds = seedpoints.size(); itk::Index<3> zeroIndex; zeroIndex.Fill(0); int 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 || stop) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { progress += print_interval; std::cout << " \r"; if (m_MaxNumTracts>0) std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << "/" << m_MaxNumTracts << '\r'; else std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << '\r'; cout.flush(); } itk::Point worldPos = 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 (IsValidPosition(worldPos)) 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; } // 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); } } counter = fib.size(); #pragma omp critical if (tractLength>=m_MinTractLength && counter>=2) { if (!stop) { - m_Tractogram.push_back(fib); + if (!m_UseOutputProbabilityMap) + m_Tractogram.push_back(fib); + else + FiberToProbmap(&fib); current_tracts++; } if (m_MaxNumTracts > 0 && current_tracts>=static_cast(m_MaxNumTracts)) { if (!stop) { std::cout << " \r"; MITK_INFO << "Reconstructed maximum number of tracts (" << current_tracts << "). Stopping tractography."; } stop = true; } } } } this->AfterTracking(); } void StreamlineTrackingFilter::CheckFiberForGmEnding(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) { ItkUcharImgType::IndexType idx; m_TissueImage->TransformPhysicalPointToIndex(fib->back(), idx); if (m_TissueImage->GetPixel(idx)==m_WmLabel) in_wm = true; else fib->pop_back(); } if (fib->size()<3 || !in_wm) { fib->clear(); return; } // 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(); } +void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib) +{ + ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0); + for (FiberType::iterator it = fib->begin(); it != fib->end(); ++it) + { + ItkDoubleImgType::IndexType idx; + itk::Point p = (*it); + m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); + + if (idx != last_idx) + { +#pragma omp critical + 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 (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() { std::cout << " \r"; - MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; - MITK_INFO << "Generating polydata "; - BuildFibers(false); + if (!m_UseOutputProbabilityMap) + { + MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; + MITK_INFO << "Generating polydata "; + BuildFibers(false); + } 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"; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h index 9d4e9ba28f..c90f06dcca 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,199 +1,206 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkMLBSTrackingFilter_h_ #define __itkMLBSTrackingFilter_h_ #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 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; } + itkGetMacro( OutputProbabilityMap, ItkDoubleImgType::Pointer) ///< Output probability map itkGetMacro( FiberPolyData, PolyDataType ) ///< Output fibers - itkSetMacro( SeedImage, ItkUcharImgType::Pointer) ///< Seeds are only placed inside of this mask. + 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( 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( 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, unsigned 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. void SetTrackingHandler( mitk::TrackingDataHandler* h ) ///< { m_TrackingHandler = h; } virtual void Update() override{ this->GenerateData(); } protected: void GenerateData() override; StreamlineTrackingFilter(); ~StreamlineTrackingFilter() {} void InitGrayMatterEndings(); - void CheckFiberForGmEnding(FiberType* fib); + void CheckFiberForGmEnding(FiberType* fib); + void FiberToProbmap(FiberType* fib); void CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir); ///< Calculate next integration step. float FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front); ///< Start streamline in one direction. bool IsValidPosition(itk::Point& pos); ///< Are we outside of the mask image? 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. float GetRandDouble(float min=-1, float max=1); 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; 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_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; ItkUcharImgType::Pointer m_StoppingRegions; ItkUcharImgType::Pointer m_SeedImage; ItkUcharImgType::Pointer m_MaskImage; ItkUcharImgType::Pointer m_TissueImage; + ItkDoubleImgType::Pointer m_OutputProbabilityMap; bool m_Verbose; bool m_AposterioriCurvCheck; bool m_AvoidStop; bool m_DemoMode; bool m_Random; + bool m_UseOutputProbabilityMap; 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; private: }; } //#ifndef ITK_MANUAL_INSTANTIATION //#include "itkMLBSTrackingFilter.cpp" //#endif #endif //__itkMLBSTrackingFilter_h_ diff --git a/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp b/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp index df1bccdb5f..5da594fb9b 100755 --- a/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp +++ b/Modules/DiffusionImaging/cmdapps/StreamlineTracking.cpp @@ -1,404 +1,431 @@ /*=================================================================== 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 #define _USE_MATH_DEFINES #include using namespace std; 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 'Tensor' algorithm)", us::Any(), false); parser.addArgument("algorithm", "a", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (Peaks, Tensor, DetODF, ProbODF, DetRF, ProbRF)", us::Any(), false); - parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle", 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("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("cutoff", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 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("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 ODF based tractography parser.addArgument("no_odf_normalization", "", mitkCommandLineParser::Bool, "Don't min-max normalize ODFs:", "No min-max normalization of ODFs"); // 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); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["input"]); string outFile = us::any_cast(parsedArgs["out"]); string algorithm = us::any_cast(parsedArgs["algorithm"]); // mitk::LoggingBackend::SetLogFile( (outFile + ".log").c_str() ); // MITK_INFO << "LOG"; 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"]); + 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"]); string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); string seedFile = ""; if (parsedArgs.count("seed_mask")) seedFile = us::any_cast(parsedArgs["seed_mask"]); string stopFile = ""; if (parsedArgs.count("stop_mask")) stopFile = us::any_cast(parsedArgs["stop_mask"]); 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 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"]); bool normalize_odfs = true; if (parsedArgs.count("no_odf_normalization")) normalize_odfs = !us::any_cast(parsedArgs["no_odf_normalization"]); 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; MITK_INFO << "loading input"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(mitk::IOUtil::LoadDataNode(input_files.at(i))->GetData()); input_images.push_back(mitkImage); } ItkUcharImgType::Pointer mask; if (!maskFile.empty()) { MITK_INFO << "loading mask image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(maskFile).GetPointer()); mask = ItkUcharImgType::New(); mitk::CastToItkImage(img, mask); } ItkUcharImgType::Pointer seed; if (!seedFile.empty()) { MITK_INFO << "loading seed image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(seedFile).GetPointer()); seed = ItkUcharImgType::New(); mitk::CastToItkImage(img, seed); } ItkUcharImgType::Pointer stop; if (!stopFile.empty()) { MITK_INFO << "loading stop image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(stopFile).GetPointer()); stop = ItkUcharImgType::New(); mitk::CastToItkImage(img, stop); } ItkUcharImgType::Pointer tissue; if (!tissueFile.empty()) { MITK_INFO << "loading tissue image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(tissueFile).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::LoadImage(file).GetPointer()); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); mitk::TrackingDataHandler* handler; if (algorithm == "DetRF" || algorithm == "ProbRF") { if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->LoadForest(forestFile); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->LoadForest(forestFile); 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 == "Tensor") { 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::Pointer 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.size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (algorithm == "DetODF" || algorithm == "ProbODF") { handler = new mitk::TrackingHandlerOdf(); for (auto input_image : input_images) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_image); caster->Update(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetOdfImage(itkImg); } dynamic_cast(handler)->SetGfaThreshold(cutoff); if (algorithm == "ProbODF") dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); if (addImages.size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); dynamic_cast(handler)->SetMinMaxNormalize(normalize_odfs); } 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->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->Update(); - vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); - mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); - if (compress>0) - outFib->Compress(compress); + 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(outFib, outFile); + mitk::IOUtil::SaveBaseData(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 363710f564..455762f504 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,406 +1,441 @@ /*=================================================================== 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 // VTK #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; QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_Controls(nullptr) { } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { } 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_StopImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TissueImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = 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_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->SetZeroEntryText("--"); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_StopImageBox->SetZeroEntryText("--"); m_Controls->m_TissueImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_TissueImageBox->SetZeroEntryText("--"); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_TissueImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); + + m_FirstTensorProbRun = true; } UpdateGui(); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer part, const QList& nodes ) { m_InputImageNodes.clear(); m_InputImages.clear(); 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())); } else if ( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); } 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())); } } } } } UpdateGui(); } void QmitkStreamlineTrackingView::UpdateGui() { m_Controls->m_TensorImageLabel->setText("mandatory"); m_Controls->m_fBox->setVisible(false); m_Controls->m_fLabel->setVisible(false); m_Controls->m_gBox->setVisible(false); m_Controls->m_gLabel->setVisible(false); m_Controls->m_FaImageBox->setVisible(false); m_Controls->mFaImageLabel->setVisible(false); m_Controls->m_NormalizeODFsBox->setVisible(false); if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_SeedGmBox->setVisible(true); else m_Controls->m_SeedGmBox->setVisible(false); if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.size()+" images selected"); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->m_InputData->setTitle("Input Data"); m_Controls->commandLinkButton->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setVisible(true); m_Controls->m_fLabel->setVisible(true); m_Controls->m_gBox->setVisible(true); m_Controls->m_gLabel->setVisible(true); m_Controls->mFaImageLabel->setVisible(true); m_Controls->m_FaImageBox->setVisible(true); if (m_Controls->m_ModeBox->currentIndex()==1) m_Controls->m_NormalizeODFsBox->setVisible(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->mFaImageLabel->setVisible(true); m_Controls->m_FaImageBox->setVisible(true); m_Controls->m_NormalizeODFsBox->setVisible(true); } else { } } else { m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->commandLinkButton->setEnabled(false); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_InputImages.empty()) return; mitk::TrackingDataHandler* trackingHandler; if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; 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!"); return; } - QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. Please keep the state of the ODF normalization box (see advanced parameters) in mind. TEND parameters are ignored."); + if (m_FirstTensorProbRun) + { + QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. Please keep the state of the ODF normalization box (see advanced parameters) in mind. TEND parameters are ignored."); + m_FirstTensorProbRun = false; + } TensorImageType::Pointer itkImg = TensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); typedef itk::TensorImageToQBallImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkImg ); filter->Update(); typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; trackingHandler = new mitk::TrackingHandlerOdf(); dynamic_cast(trackingHandler)->SetOdfImage(filter->GetOutput()); dynamic_cast(trackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(trackingHandler)->SetMinMaxNormalize(m_Controls->m_NormalizeODFsBox->isChecked()); dynamic_cast(trackingHandler)->SetOdfPower(1); 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(trackingHandler)->SetGfaImage(itkImg); } } else { trackingHandler = new mitk::TrackingHandlerTensor(); for (int i=0; i<(int)m_InputImages.size(); i++) { TensorImageType::Pointer itkImg = TensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(i), itkImg); dynamic_cast(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(trackingHandler)->SetFaImage(itkImg); } dynamic_cast(trackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(trackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(trackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; trackingHandler = new mitk::TrackingHandlerOdf(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); dynamic_cast(trackingHandler)->SetOdfImage(itkImg); dynamic_cast(trackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(trackingHandler)->SetMinMaxNormalize(m_Controls->m_NormalizeODFsBox->isChecked()); 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(trackingHandler)->SetGfaImage(itkImg); } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is only implementedfor ODF images."); return; } try { 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(); trackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(trackingHandler)->SetPeakImage(itkImg); dynamic_cast(trackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); } catch(...) { return; } } trackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); trackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); trackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); trackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: trackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if (m_Controls->m_SeedImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_SeedImageBox->GetSelectedNode()->GetData()), mask); tracker->SetSeedImage(mask); } if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); tracker->SetMaskImage(mask); } if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); tracker->SetStoppingRegions(mask); } if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TissueImageBox->GetSelectedNode()->GetData()), mask); tracker->SetTissueImage(mask); tracker->SetSeedOnlyGm(m_Controls->m_SeedGmBox->isChecked()); } tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); //tracker->SetSamplingDistance(0.7); tracker->SetUseStopVotes(false); tracker->SetOnlyForwardSamples(false); tracker->SetAposterioriCurvCheck(false); tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); tracker->SetTrackingHandler(trackingHandler); tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); - tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); + tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); + tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); tracker->Update(); delete trackingHandler; - vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); - if ( fiberBundle->GetNumberOfLines()==0 ) - { - QMessageBox warnBox; - warnBox.setWindowTitle("Warning"); - warnBox.setText("No fiberbundle was generated!"); - warnBox.setDetailedText("No fibers were generated using the parameters: \n\n" + m_Controls->m_FaThresholdLabel->text() + "\n" + m_Controls->m_AngularThresholdLabel->text() + "\n" + m_Controls->m_fLabel->text() + "\n" + m_Controls->m_gLabel->text() + "\n" + m_Controls->m_StepsizeLabel->text() + "\n" + m_Controls->m_MinTractLengthLabel->text() + "\n" + m_Controls->m_SeedsPerVoxelLabel->text() + "\n\nPlease check your parametersettings."); - warnBox.setIcon(QMessageBox::Warning); - warnBox.exec(); - return; - } - mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); - fib->SetReferenceGeometry(dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()); - if (m_Controls->m_ResampleFibersBox->isChecked()) - fib->Compress(m_Controls->m_FiberErrorBox->value()); - fib->ColorFibersByOrientation(); - - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(fib); - QString name("FiberBundle_"); - name += m_InputImageNodes.at(0)->GetName().c_str(); - name += "_Streamline"; - node->SetName(name.toStdString()); - node->SetVisibility(true); - GetDataStorage()->Add(node, m_InputImageNodes.at(0)); + if (!tracker->GetUseOutputProbabilityMap()) + { + vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); + if (fiberBundle->GetNumberOfLines() == 0) + { + QMessageBox warnBox; + warnBox.setWindowTitle("Warning"); + warnBox.setText("No fiberbundle was generated!"); + warnBox.setDetailedText("No fibers were generated using the parameters: \n\n" + m_Controls->m_FaThresholdLabel->text() + "\n" + m_Controls->m_AngularThresholdLabel->text() + "\n" + m_Controls->m_fLabel->text() + "\n" + m_Controls->m_gLabel->text() + "\n" + m_Controls->m_StepsizeLabel->text() + "\n" + m_Controls->m_MinTractLengthLabel->text() + "\n" + m_Controls->m_SeedsPerVoxelLabel->text() + "\n\nPlease check your parametersettings."); + warnBox.setIcon(QMessageBox::Warning); + warnBox.exec(); + return; + } + mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); + fib->SetReferenceGeometry(dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()); + if (m_Controls->m_ResampleFibersBox->isChecked()) + fib->Compress(m_Controls->m_FiberErrorBox->value()); + fib->ColorFibersByOrientation(); + + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(fib); + QString name("FiberBundle_"); + name += m_InputImageNodes.at(0)->GetName().c_str(); + name += "_Streamline"; + node->SetName(name.toStdString()); + node->SetVisibility(true); + GetDataStorage()->Add(node, m_InputImageNodes.at(0)); + } + else + { + TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + QString name("ProbabilityMap_"); + name += m_InputImageNodes.at(0)->GetName().c_str(); + node->SetName(name.toStdString()); + node->SetVisibility(true); + + 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_InputImageNodes.at(0)); + } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h index 7c228ed9fd..ac09642ea2 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.h @@ -1,84 +1,85 @@ /*=================================================================== 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 QmitkStreamlineTrackingView_h #define QmitkStreamlineTrackingView_h #include #include "ui_QmitkStreamlineTrackingViewControls.h" #include #include #include #include #include #include #include #include /*! \brief View for tensor based deterministic streamline fiber tracking. */ class QmitkStreamlineTrackingView : public QmitkAbstractView { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; typedef itk::Image< float, 3 > ItkFloatImageType; QmitkStreamlineTrackingView(); virtual ~QmitkStreamlineTrackingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; protected slots: void DoFiberTracking(); ///< start fiber tracking void UpdateGui(); protected: /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkStreamlineTrackingViewControls* m_Controls; protected slots: private: std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input images std::vector< mitk::Image::Pointer > m_InputImages; ///< input image datanode + bool m_FirstTensorProbRun; }; #endif // _QMITKFIBERTRACKINGVIEW_H_INCLUDED 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 0c52cb6df3..d0d2631389 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,672 +1,685 @@ QmitkStreamlineTrackingViewControls 0 0 358 733 0 0 QmitkTemplate 0 0 3 3 false Start Tractography Qt::Vertical QSizePolicy::Expanding 20 220 Advanced Paramaters Step Size: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Angular Threshold: Default: 90° * step_size -1 90 1 -1 Min. Tract Length: Shorter fibers are discarded. 1 999.000000000000000 1.000000000000000 20.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: 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 g: 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 Flip directions: 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 Neighborhood Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 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 - - + + - + If false, nearest neighbor interpolation is used. - Min-max normalize ODFs + Enable Trilinear Interpolation true - - + + - If false, nearest neighbor interpolation is used. + Requires tissue image. - Enable Trilinear Interpolation + Enable gray matter seeding + + + false + + + + + + + + + + Min-max normalize ODFs true - - + + - Requires tissue image. + Output probability map instead of tractogram. - Enable gray matter seeding + Output probability map false Please Select Input Data Only track insida mask area. Mask Image: Fibers that enter a region defined in this image will stop immediately. QComboBox::AdjustToMinimumContentsLength - If an image is selected, the stopping criterion is not calculated from the input image but instead the selected image is used. QComboBox::AdjustToMinimumContentsLength - Binary seed ROI. If not specified, the whole image area is seeded. Seed Image: <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Seed points are only placed inside the mask image. QComboBox::AdjustToMinimumContentsLength - Only track insida mask area. Stop Image: Input DTI - Tensor/Peak Image: + Tensor/Peak/ODF Image: Only track insida mask area. FA/GFA Image: Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. QComboBox::AdjustToMinimumContentsLength - Only track insida mask area. Tissue Image: Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. QComboBox::AdjustToMinimumContentsLength - 0 0 Parameters 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. 999999999 Cutoff: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 Qt::Horizontal QSizePolicy::Fixed 200 0 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h