diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 6f42d82d64..7ef5f29891 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,994 +1,993 @@ /*=================================================================== 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) + , m_UseOutputProbabilityMap(false) { 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; } if (m_SeedPoints.empty()) GetSeedPointsFromSeedImage(); m_BuildFibersReady = 0; m_BuildFibersFinished = false; m_Tractogram.clear(); m_SamplingPointset = mitk::PointSet::New(); m_AlternativePointset = mitk::PointSet::New(); m_StopVotePointset = mitk::PointSet::New(); m_StartTime = std::chrono::system_clock::now(); if (m_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; } bool StreamlineTrackingFilter ::IsInGm(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; } 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 && !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::GetSeedPointsFromSeedImage() { MITK_INFO << "Calculating seed points."; m_SeedPoints.clear(); if (!m_ControlGmEndings) { typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); sit.GoToBegin(); while (!sit.IsAtEnd()) { if (sit.Value()>0) { 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 (IsValidPosition(worldPos) && (!m_SeedOnlyGm || IsInGm(worldPos))) { m_SeedPoints.push_back(worldPos); for (int s = 1; s < m_SeedsPerVoxel; s++) { 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); 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]); } } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (m_Random) { - std::srand(std::time(0)); - std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); + std::srand(std::time(0)); + std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); } bool stop = false; unsigned int current_tracts = 0; - int num_seeds = m_SeedPoints.size(); + int num_seeds = m_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 = m_SeedPoints.at(temp_i); + 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 (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) { - if (!m_UseOutputProbabilityMap) - m_Tractogram.push_back(fib); - else - FiberToProbmap(&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 last_idx; last_idx.Fill(0); + for (auto p : *fib) { - ItkDoubleImgType::IndexType idx; - itk::Point p = (*it); + ItkDoubleImgType::IndexType idx; m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); - if (idx != last_idx) - { -#pragma omp critical - if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(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 (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"; + if (m_Verbose) + std::cout << " \r"; 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"; + 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(); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp index 006c3a89f0..8135c7dcf4 100755 --- a/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Testing/mitkStreamlineTrackingTest.cpp @@ -1,138 +1,137 @@ /*=================================================================== 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 using namespace std; int mitkStreamlineTrackingTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkStreamlineTrackingTest"); MITK_TEST_CONDITION_REQUIRED(argc>3,"check for input data") string dtiFileName = argv[1]; string maskFileName = argv[2]; string referenceFileName = argv[3]; MITK_INFO << "DTI file: " << dtiFileName; MITK_INFO << "Mask file: " << maskFileName; MITK_INFO << "Reference fiber file: " << referenceFileName; float minFA = 0.05; - float minCurv = -1; float stepSize = -1; float tendf = 1; float tendg = 0; float minLength = 20; int numSeeds = 1; bool interpolate = true; try { MITK_INFO << "Loading tensor image ..."; typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; mitk::TensorImage::Pointer mitkTensorImage = dynamic_cast(mitk::IOUtil::LoadDataNode(dtiFileName)->GetData()); ItkTensorImage::Pointer itk_dti = ItkTensorImage::New(); mitk::CastToItkImage(mitkTensorImage, itk_dti); MITK_INFO << "Loading seed image ..."; typedef itk::Image< unsigned char, 3 > ItkUCharImageType; mitk::Image::Pointer mitkSeedImage = mitk::IOUtil::LoadImage(maskFileName); MITK_INFO << "Loading mask image ..."; mitk::Image::Pointer mitkMaskImage = mitk::IOUtil::LoadImage(maskFileName); omp_set_num_threads(1); // tensor tracking handler mitk::TrackingHandlerTensor* handler = new mitk::TrackingHandlerTensor(); handler->SetF(tendf); handler->SetG(tendg); handler->SetInterpolate(interpolate); handler->SetFaThreshold(minFA); handler->AddTensorImage(itk_dti); // tracker typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); tracker->SetSeedsPerVoxel(numSeeds); tracker->SetNumberOfSamples(0); tracker->SetStepSize(stepSize); tracker->SetAposterioriCurvCheck(false); tracker->SetSeedOnlyGm(false); tracker->SetTrackingHandler(handler); tracker->SetMinTractLength(minLength); tracker->SetRandom(false); if (mitkSeedImage.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkSeedImage, mask); tracker->SetSeedImage(mask); } if (mitkMaskImage.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(mitkMaskImage, mask); tracker->SetMaskImage(mask); } tracker->Update(); vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer fib1 = mitk::FiberBundle::New(fiberBundle); mitk::FiberBundle::Pointer fib2 = dynamic_cast(mitk::IOUtil::LoadDataNode(referenceFileName)->GetData()); MITK_TEST_CONDITION_REQUIRED(fib2.IsNotNull(), "Check if reference tractogram is not null."); bool ok = fib1->Equals(fib2); if (!ok) { MITK_WARN << "TEST FAILED. TRACTOGRAMS ARE NOT EQUAL!"; mitk::IOUtil::SaveBaseData(fib1, mitk::IOUtil::GetTempPath()+"testBundle.fib"); MITK_INFO << "OUTPUT: " << mitk::IOUtil::GetTempPath(); } MITK_TEST_CONDITION_REQUIRED(ok, "Check if tractograms are equal."); delete handler; } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_TEST_END(); }