diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp index d1d996a845..111bf7baf6 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp @@ -1,285 +1,287 @@ /*=================================================================== 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 "mitkTrackingHandlerPeaks.h" namespace mitk { TrackingHandlerPeaks::TrackingHandlerPeaks() : m_PeakThreshold(0.1) , m_ApplyDirectionMatrix(false) { } TrackingHandlerPeaks::~TrackingHandlerPeaks() { } void TrackingHandlerPeaks::InitForTracking() { MITK_INFO << "Initializing peak tracker."; if (m_NeedsDataInit) { itk::Vector spacing4 = m_PeakImage->GetSpacing(); itk::Point origin4 = m_PeakImage->GetOrigin(); itk::Matrix direction4 = m_PeakImage->GetDirection(); itk::ImageRegion<4> imageRegion4 = m_PeakImage->GetLargestPossibleRegion(); spacing3[0] = spacing4[0]; spacing3[1] = spacing4[1]; spacing3[2] = spacing4[2]; origin3[0] = origin4[0]; origin3[1] = origin4[1]; origin3[2] = origin4[2]; for (int r=0; r<3; r++) for (int c=0; c<3; c++) { direction3[r][c] = direction4[r][c]; m_FloatImageRotation[r][c] = direction4[r][c]; } imageRegion3.SetSize(0, imageRegion4.GetSize()[0]); imageRegion3.SetSize(1, imageRegion4.GetSize()[1]); imageRegion3.SetSize(2, imageRegion4.GetSize()[2]); m_DummyImage = ItkUcharImgType::New(); m_DummyImage->SetSpacing( spacing3 ); m_DummyImage->SetOrigin( origin3 ); m_DummyImage->SetDirection( direction3 ); m_DummyImage->SetRegions( imageRegion3 ); m_DummyImage->Allocate(); m_DummyImage->FillBuffer(0.0); m_NumDirs = imageRegion4.GetSize(3)/3; m_NeedsDataInit = false; } std::cout << "TrackingHandlerPeaks - Peak threshold: " << m_PeakThreshold << std::endl; } vnl_vector_fixed TrackingHandlerPeaks::GetMatchingDirection(itk::Index<3> idx3, vnl_vector_fixed& oldDir) { vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0; float mag = oldDir.magnitude(); if (magGetIntegerVariate(m_NumDirs-1); + int i = 0; +#pragma omp critical + i = m_RngItk->GetIntegerVariate(m_NumDirs-1); out_dir = GetDirection(idx3, i); if (out_dir.magnitude()>mitk::eps) { oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; found = true; break; } } } if (!found) { // if you didn't find a non-zero random direction, take first non-zero direction you find for (int i=0; imitk::eps) { oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } } } else { for (int i=0; i dir = GetDirection(idx3, i); mag = dir.magnitude(); if (mag>mitk::eps) dir.normalize(); float a = dot_product(dir, oldDir); if (fabs(a)>angle) { angle = fabs(a); if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= mag; out_dir *= angle; // shrink contribution of direction if is less parallel to previous direction } } } return out_dir; } vnl_vector_fixed TrackingHandlerPeaks::GetDirection(itk::Index<3> idx3, int dirIdx) { vnl_vector_fixed dir; dir.fill(0.0); if ( !m_DummyImage->GetLargestPossibleRegion().IsInside(idx3) ) return dir; PeakImgType::IndexType idx4; idx4.SetElement(0,idx3[0]); idx4.SetElement(1,idx3[1]); idx4.SetElement(2,idx3[2]); for (int k=0; k<3; k++) { idx4.SetElement(3, dirIdx*3 + k); dir[k] = m_PeakImage->GetPixel(idx4); } if (m_FlipX) dir[0] *= -1; if (m_FlipY) dir[1] *= -1; if (m_FlipZ) dir[2] *= -1; if (m_ApplyDirectionMatrix) dir = m_FloatImageRotation*dir; return dir; } vnl_vector_fixed TrackingHandlerPeaks::GetDirection(itk::Point itkP, bool interpolate, vnl_vector_fixed oldDir){ // transform physical point to index coordinates itk::Index<3> idx3; itk::ContinuousIndex< float, 3> cIdx; m_DummyImage->TransformPhysicalPointToIndex(itkP, idx3); m_DummyImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); vnl_vector_fixed dir; dir.fill(0.0); if ( !m_DummyImage->GetLargestPossibleRegion().IsInside(idx3) ) return dir; if (interpolate) { float frac_x = cIdx[0] - idx3[0]; float frac_y = cIdx[1] - idx3[1]; float frac_z = cIdx[2] - idx3[2]; if (frac_x<0) { idx3[0] -= 1; frac_x += 1; } if (frac_y<0) { idx3[1] -= 1; frac_y += 1; } if (frac_z<0) { idx3[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (idx3[0] >= 0 && idx3[0] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(0) - 1) && idx3[1] >= 0 && idx3[1] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(1) - 1) && idx3[2] >= 0 && idx3[2] < static_cast(m_DummyImage->GetLargestPossibleRegion().GetSize(2) - 1)) { // trilinear interpolation vnl_vector_fixed interpWeights; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); dir = GetMatchingDirection(idx3, oldDir) * interpWeights[0]; itk::Index<3> tmpIdx = idx3; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[1]; tmpIdx = idx3; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[2]; tmpIdx = idx3; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[3]; tmpIdx = idx3; tmpIdx[0]++; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[4]; tmpIdx = idx3; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[5]; tmpIdx = idx3; tmpIdx[2]++; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[6]; tmpIdx = idx3; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir) * interpWeights[7]; } } else dir = GetMatchingDirection(idx3, oldDir); return dir; } vnl_vector_fixed TrackingHandlerPeaks::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { // CHECK: wann wird wo normalisiert vnl_vector_fixed output_direction; output_direction.fill(0); itk::Index<3> index; m_DummyImage->TransformPhysicalPointToIndex(pos, index); vnl_vector_fixed oldDir = olddirs.back(); float old_mag = oldDir.magnitude(); if (!m_Interpolate && oldIndex==index) return oldDir; output_direction = GetDirection(pos, m_Interpolate, oldDir); float mag = output_direction.magnitude(); if (mag>=m_PeakThreshold) { output_direction.normalize(); float a = 1; if (old_mag>0.5) a = dot_product(output_direction, oldDir); if (a>=m_AngularThreshold) output_direction *= mag; else output_direction.fill(0); } else output_direction.fill(0); return output_direction; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 942bc31431..0350464ee2 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1004 +1,1020 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include "itkStreamlineTrackingFilter.h" #include #include #include #include #include "itkPointShell.h" #include +#include #define _USE_MATH_DEFINES #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_BuildFibersFinished(false) , m_BuildFibersReady(0) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_StoppingRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_TissueImage(nullptr) , m_OutputProbabilityMap(nullptr) , m_AngularThresholdDeg(-1) , m_StepSizeVox(-1) , m_SamplingDistanceVox(-1) , m_AngularThreshold(-1) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) , m_AvoidStop(true) , m_RandomSampling(false) , m_SamplingDistance(-1) , m_DeflectionMod(1.0) , m_OnlyForwardSamples(true) , m_UseStopVotes(true) , m_NumberOfSamples(30) , m_NumPreviousDirections(1) , m_WmLabel(3) // mrtrix 5ttseg labels , m_GmLabel(1) // mrtrix 5ttseg labels , m_SeedOnlyGm(false) , m_ControlGmEndings(false) , m_MaxNumTracts(-1) , m_Verbose(true) , m_AposterioriCurvCheck(false) , m_DemoMode(false) , m_Random(true) , m_UseOutputProbabilityMap(false) + , m_CurrentTracts(0) + , m_Progress(0) + , m_StopTracking(false) { this->SetNumberOfRequiredInputs(0); } +std::string StreamlineTrackingFilter::GetStatusText() +{ + std::string status = "Seedpoints processed: " + boost::lexical_cast(m_Progress) + "/" + boost::lexical_cast(m_SeedPoints.size()); + if (m_SeedPoints.size()>0) + status += " (" + boost::lexical_cast(100*m_Progress/m_SeedPoints.size()) + "%)"; + if (m_MaxNumTracts>0) + status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts) + "/" + boost::lexical_cast(m_MaxNumTracts); + else + status += "\nFibers accepted: " + boost::lexical_cast(m_CurrentTracts); + + return status; +} void StreamlineTrackingFilter::BeforeTracking() { + m_StopTracking = false; m_TrackingHandler->SetRandom(m_Random); m_TrackingHandler->InitForTracking(); m_FiberPolyData = PolyDataType::New(); m_Points = vtkSmartPointer< vtkPoints >::New(); m_Cells = vtkSmartPointer< vtkCellArray >::New(); itk::Vector< double, 3 > imageSpacing = m_TrackingHandler->GetSpacing(); float minSpacing; if(imageSpacing[0]SetAngularThreshold(m_AngularThreshold); if (m_SamplingDistanceVoxGetNumberOfThreads(); i++) { PolyDataType poly = PolyDataType::New(); m_PolyDataContainer.push_back(poly); } if (m_UseOutputProbabilityMap) { m_OutputProbabilityMap = ItkDoubleImgType::New(); m_OutputProbabilityMap->SetSpacing(imageSpacing); m_OutputProbabilityMap->SetOrigin(m_TrackingHandler->GetOrigin()); m_OutputProbabilityMap->SetDirection(m_TrackingHandler->GetDirection()); m_OutputProbabilityMap->SetRegions(m_TrackingHandler->GetLargestPossibleRegion()); m_OutputProbabilityMap->Allocate(); m_OutputProbabilityMap->FillBuffer(0); } 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_SeedPoints.empty()) GetSeedPointsFromSeedImage(); else m_SeedOnlyGm = false; if (m_TissueImage.IsNull()) { if (m_SeedOnlyGm) { MITK_WARN << "StreamlineTracking - Cannot seed in gray matter. No tissue type image set."; m_SeedOnlyGm = false; } if (m_ControlGmEndings) { MITK_WARN << "StreamlineTracking - Cannot control gray matter endings. No tissue type image set."; m_ControlGmEndings = false; } } else { if (m_ControlGmEndings) m_SeedOnlyGm = true; if (m_ControlGmEndings || m_SeedOnlyGm) std::cout << "StreamlineTracking - Using tissue image" << std::endl; else MITK_WARN << "StreamlineTracking - Tissue image set but no gray matter seeding or fiber endpoint-control enabled" << std::endl; } m_BuildFibersReady = 0; m_BuildFibersFinished = false; m_Tractogram.clear(); m_SamplingPointset = mitk::PointSet::New(); m_AlternativePointset = mitk::PointSet::New(); m_StopVotePointset = mitk::PointSet::New(); m_StartTime = std::chrono::system_clock::now(); if (m_SeedOnlyGm && m_ControlGmEndings) InitGrayMatterEndings(); if (m_DemoMode) omp_set_num_threads(1); if (m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::DETERMINISTIC) std::cout << "StreamlineTracking - Mode: deterministic" << std::endl; else if(m_TrackingHandler->GetMode()==mitk::TrackingDataHandler::MODE::PROBABILISTIC) std::cout << "StreamlineTracking - Mode: probabilistic" << std::endl; else std::cout << "StreamlineTracking - Mode: ???" << std::endl; std::cout << "StreamlineTracking - Angular threshold: " << m_AngularThreshold << " (" << 180*std::acos( m_AngularThreshold )/M_PI << "°)" << std::endl; std::cout << "StreamlineTracking - Stepsize: " << m_StepSize << "mm (" << m_StepSize/minSpacing << "*vox)" << std::endl; std::cout << "StreamlineTracking - Seeds per voxel: " << m_SeedsPerVoxel << std::endl; std::cout << "StreamlineTracking - Max. tract length: " << m_MaxTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Min. tract length: " << m_MinTractLength << "mm" << std::endl; std::cout << "StreamlineTracking - Max. num. tracts: " << m_MaxNumTracts << std::endl; std::cout << "StreamlineTracking - Num. neighborhood samples: " << m_NumberOfSamples << std::endl; std::cout << "StreamlineTracking - Max. sampling distance: " << m_SamplingDistance << "mm (" << m_SamplingDistance/minSpacing << "*vox)" << std::endl; std::cout << "StreamlineTracking - Deflection modifier: " << m_DeflectionMod << std::endl; std::cout << "StreamlineTracking - Use stop votes: " << m_UseStopVotes << std::endl; std::cout << "StreamlineTracking - Only frontal samples: " << m_OnlyForwardSamples << std::endl; if (m_DemoMode) { std::cout << "StreamlineTracking - Running in demo mode"; std::cout << "StreamlineTracking - Starting streamline tracking using 1 thread" << std::endl; } else std::cout << "StreamlineTracking - Starting streamline tracking using " << omp_get_max_threads() << " threads" << std::endl; } void StreamlineTrackingFilter::InitGrayMatterEndings() { m_TrackingHandler->SetAngularThreshold(0); m_GmStubs.clear(); if (m_TissueImage.IsNotNull()) { std::cout << "StreamlineTracking - initializing GM endings" << std::endl; ImageRegionConstIterator< ItkUcharImgType > it(m_TissueImage, m_TissueImage->GetLargestPossibleRegion() ); it.GoToBegin(); vnl_vector_fixed d1; d1.fill(0.0); std::deque< vnl_vector_fixed > olddirs; while (olddirs.size() start; m_TissueImage->TransformIndexToPhysicalPoint(s_idx, start); itk::Point wm_p; float max = -1; FiberType fib; for (int x : {-1,0,1}) for (int y : {-1,0,1}) for (int z : {-1,0,1}) { if (x==y && y==z) continue; ItkUcharImgType::IndexType e_idx; e_idx[0] = s_idx[0] + x; e_idx[1] = s_idx[1] + y; e_idx[2] = s_idx[2] + z; if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_WmLabel ) continue; itk::ContinuousIndex end; m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); d1 = m_TrackingHandler->ProposeDirection(end, olddirs, s_idx); if (d1.magnitude()<0.0001) continue; d1.normalize(); vnl_vector_fixed< float, 3 > d2; d2[0] = end[0] - start[0]; d2[1] = end[1] - start[1]; d2[2] = end[2] - start[2]; d2.normalize(); float a = fabs(dot_product(d1,d2)); if (a>max) { max = a; wm_p = end; } } if (max>=0) { fib.push_back(start); fib.push_back(wm_p); m_GmStubs.push_back(fib); } } ++it; } } m_TrackingHandler->SetAngularThreshold(m_AngularThreshold); } void StreamlineTrackingFilter::CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir) { pos[0] += dir[0]*m_StepSize; pos[1] += dir[1]*m_StepSize; pos[2] += dir[2]*m_StepSize; } bool StreamlineTrackingFilter ::IsValidPosition(const 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(const itk::Point &pos) { if (m_TissueImage.IsNull()) return true; ItkUcharImgType::IndexType idx; m_TissueImage->TransformPhysicalPointToIndex(pos, idx); if (m_TissueImage->GetLargestPossibleRegion().IsInside(idx) && m_TissueImage->GetPixel(idx) == m_GmLabel) return true; return false; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< float > theta; theta.resize(NPoints); std::vector< float > phi; phi.resize(NPoints); float C = sqrt(4*M_PI); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(int i=0; i0 && i d; d[0] = cos(theta[i]) * cos(phi[i]); d[1] = cos(theta[i]) * sin(phi[i]); d[2] = sin(theta[i]); pointshell.push_back(d); } return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { if (m_DemoMode) { m_SamplingPointset->Clear(); m_AlternativePointset->Clear(); m_StopVotePointset->Clear(); } vnl_vector_fixed direction; direction.fill(0); 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 (unsigned int i=0; i d; bool is_stop_voter = false; if (m_Random && m_RandomSampling) { d[0] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[1] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d[2] = m_TrackingHandler->GetRandDouble(-0.5, 0.5); d.normalize(); d *= m_TrackingHandler->GetRandDouble(0,m_SamplingDistance); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); if (m_UseStopVotes && dot>0.7) { is_stop_voter = true; possible_stop_votes++; } else if (m_OnlyForwardSamples && dot<0) continue; d *= m_SamplingDistance; } sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (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; } if (m_DemoMode && !m_UseOutputProbabilityMap) // CHECK: warum sind die samplingpunkte der streamline in der visualisierung immer einen schritt voras? { #pragma omp critical { m_BuildFibersReady++; m_Tractogram.push_back(*fib); BuildFibers(true); m_Stop = true; while (m_Stop){ } } } dir.normalize(); last_dirs.push_back(dir); if (last_dirs.size()>m_NumPreviousDirections) last_dirs.pop_front(); dir = GetNewDirection(pos, last_dirs, oldIndex); while (m_PauseTracking){} if (dir.magnitude()<0.0001) return tractLength; } return tractLength; } int StreamlineTrackingFilter::CheckCurvature(FiberType* fib, bool front) { float m_Distance = 5; if (fib->size()<3) return 0; float dist = 0; std::vector< vnl_vector_fixed< float, 3 > > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0); float dev = 0; if (front) { int c=0; while(distsize()-1) { itk::Point p1 = fib->at(c); itk::Point p2 = fib->at(c+1); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==0) meanV += v; c++; } } else { int c=fib->size()-1; while(dist0) { itk::Point p1 = fib->at(c); itk::Point p2 = fib->at(c-1); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==(int)fib->size()-1) meanV += v; c--; } } meanV.normalize(); for (unsigned int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); if (dev<30) return 0; else return vectors.size(); } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { MITK_INFO << "StreamlineTracking - Calculating seed points."; m_SeedPoints.clear(); if (!m_ControlGmEndings) { typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; MaskIteratorType sit(m_SeedImage, m_SeedImage->GetLargestPossibleRegion()); sit.GoToBegin(); 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] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[1] = index[1] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); start[2] = index[2] + m_TrackingHandler->GetRandDouble(-0.5, 0.5); itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); m_SeedPoints.push_back(worldPos); } } } ++sit; } } else { for (auto s : m_GmStubs) m_SeedPoints.push_back(s[1]); } } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (m_Random) std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); - bool stop = false; - unsigned int current_tracts = 0; + m_CurrentTracts = 0; int num_seeds = m_SeedPoints.size(); itk::Index<3> zeroIndex; zeroIndex.Fill(0); - int progress = 0; + m_Progress = 0; int i = 0; int print_interval = num_seeds/100; if (print_interval<100) m_Verbose=false; #pragma omp parallel - while (i=num_seeds || stop) + if (temp_i>=num_seeds || m_StopTracking) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { - progress += print_interval; + m_Progress += print_interval; std::cout << " \r"; if (m_MaxNumTracts>0) - std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << "/" << m_MaxNumTracts << '\r'; + std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << "/" << m_MaxNumTracts << '\r'; else - std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << '\r'; + std::cout << "Tried: " << m_Progress << "/" << num_seeds << " | Accepted: " << m_CurrentTracts << '\r'; cout.flush(); } const itk::Point worldPos = m_SeedPoints.at(temp_i); FiberType fib; float tractLength = 0; unsigned int counter = 0; // get starting direction vnl_vector_fixed dir; dir.fill(0.0); std::deque< vnl_vector_fixed > olddirs; while (olddirs.size() gm_start_dir; if (m_ControlGmEndings) { gm_start_dir[0] = m_GmStubs[temp_i][1][0] - m_GmStubs[temp_i][0][0]; gm_start_dir[1] = m_GmStubs[temp_i][1][1] - m_GmStubs[temp_i][0][1]; gm_start_dir[2] = m_GmStubs[temp_i][1][2] - m_GmStubs[temp_i][0][2]; gm_start_dir.normalize(); olddirs.pop_back(); olddirs.push_back(gm_start_dir); } if (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_StopTracking) { if (!m_UseOutputProbabilityMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); - current_tracts++; + m_CurrentTracts++; } - if (m_MaxNumTracts > 0 && current_tracts>=static_cast(m_MaxNumTracts)) + if (m_MaxNumTracts > 0 && m_CurrentTracts>=static_cast(m_MaxNumTracts)) { - if (!stop) + if (!m_StopTracking) { std::cout << " \r"; - MITK_INFO << "Reconstructed maximum number of tracts (" << current_tracts << "). Stopping tractography."; + MITK_INFO << "Reconstructed maximum number of tracts (" << m_CurrentTracts << "). Stopping tractography."; } - stop = true; + m_StopTracking = 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 (auto p : *fib) { ItkDoubleImgType::IndexType idx; m_OutputProbabilityMap->TransformPhysicalPointToIndex(p, idx); if (idx != last_idx) { if (m_OutputProbabilityMap->GetLargestPossibleRegion().IsInside(idx)) m_OutputProbabilityMap->SetPixel(idx, m_OutputProbabilityMap->GetPixel(idx)+1); last_idx = idx; } } } void StreamlineTrackingFilter::BuildFibers(bool check) { if (m_BuildFibersReady::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); for (unsigned int i=0; i container = vtkSmartPointer::New(); FiberType fib = m_Tractogram.at(i); for (FiberType::iterator it = fib.begin(); it!=fib.end(); ++it) { vtkIdType id = vNewPoints->InsertNextPoint((*it).GetDataPointer()); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if (check) for (int i=0; iSetPoints(vNewPoints); m_FiberPolyData->SetLines(vNewLines); m_BuildFibersFinished = true; } void StreamlineTrackingFilter::AfterTracking() { if (m_Verbose) std::cout << " \r"; if (!m_UseOutputProbabilityMap) { MITK_INFO << "Reconstructed " << m_Tractogram.size() << " fibers."; MITK_INFO << "Generating polydata "; BuildFibers(false); } else { itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::Pointer filter = itk::RescaleIntensityImageFilter< ItkDoubleImgType, ItkDoubleImgType >::New(); filter->SetInput(m_OutputProbabilityMap); filter->SetOutputMaximum(1.0); filter->SetOutputMinimum(0.0); filter->Update(); m_OutputProbabilityMap = filter->GetOutput(); } MITK_INFO << "done"; m_EndTime = std::chrono::system_clock::now(); std::chrono::hours hh = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::minutes mm = std::chrono::duration_cast(m_EndTime - m_StartTime); std::chrono::seconds ss = std::chrono::duration_cast(m_EndTime - m_StartTime); mm %= 60; ss %= 60; MITK_INFO << "Tracking took " << hh.count() << "h, " << mm.count() << "m and " << ss.count() << "s"; m_SeedPoints.clear(); } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h index 8a56fdca7c..7a74725311 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.h @@ -1,216 +1,222 @@ /*=================================================================== 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 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. + itkSetMacro( StopTracking, bool) ///< Set flag to stop tractography void SetSeedPoints(std::vector< itk::Point > seedPoints) ///< Use manually defined points in physical space as seed points instead of seed image { m_SeedPoints = seedPoints; this->Modified(); } void SetTrackingHandler( mitk::TrackingDataHandler* h ) ///< { m_TrackingHandler = h; } virtual void Update() override{ this->GenerateData(); } + std::string GetStatusText(); + protected: void GenerateData() override; StreamlineTrackingFilter(); ~StreamlineTrackingFilter() {} void InitGrayMatterEndings(); void CheckFiberForGmEnding(FiberType* fib); void FiberToProbmap(FiberType* fib); void GetSeedPointsFromSeedImage(); void CalculateNewPosition(itk::Point& pos, vnl_vector_fixed& dir); ///< Calculate next integration step. float FollowStreamline(itk::Point start_pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front); ///< Start streamline in one direction. bool IsValidPosition(const itk::Point& pos); ///< Are we outside of the mask image? bool IsInGm(const itk::Point &pos); vnl_vector_fixed GetNewDirection(itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< Determine new direction by sample voting at the current position taking the last progression direction into account. std::vector< vnl_vector_fixed > CreateDirections(int NPoints); void BeforeTracking(); void AfterTracking(); PolyDataType m_FiberPolyData; vtkSmartPointer m_Points; vtkSmartPointer m_Cells; BundleType m_Tractogram; BundleType m_GmStubs; ItkUcharImgType::Pointer m_StoppingRegions; ItkUcharImgType::Pointer m_SeedImage; ItkUcharImgType::Pointer m_MaskImage; ItkUcharImgType::Pointer m_TissueImage; ItkDoubleImgType::Pointer m_OutputProbabilityMap; float m_AngularThresholdDeg; float m_StepSizeVox; float m_SamplingDistanceVox; float m_AngularThreshold; float m_StepSize; int m_MaxLength; float m_MinTractLength; float m_MaxTractLength; int m_SeedsPerVoxel; bool m_AvoidStop; bool m_RandomSampling; float m_SamplingDistance; float m_DeflectionMod; bool m_OnlyForwardSamples; bool m_UseStopVotes; unsigned int m_NumberOfSamples; unsigned int m_NumPreviousDirections; int m_WmLabel; int m_GmLabel; bool m_SeedOnlyGm; bool m_ControlGmEndings; int m_MaxNumTracts; bool m_Verbose; bool m_AposterioriCurvCheck; bool m_DemoMode; bool m_Random; bool m_UseOutputProbabilityMap; std::vector< itk::Point > m_SeedPoints; + unsigned int m_CurrentTracts; + unsigned int m_Progress; + bool m_StopTracking; 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/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp index 8d7c3871fe..9418e0878b 100755 --- a/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp +++ b/Modules/DiffusionImaging/FiberTracking/IODataStructures/FiberBundle/mitkFiberBundle.cpp @@ -1,2329 +1,2328 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #define _USE_MATH_DEFINES #include "mitkFiberBundle.h" #include #include #include #include "mitkImagePixelReadAccessor.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include const char* mitk::FiberBundle::FIBER_ID_ARRAY = "Fiber_IDs"; using namespace std; mitk::FiberBundle::FiberBundle( vtkPolyData* fiberPolyData ) : m_NumFibers(0) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberPolyData = vtkSmartPointer::New(); if (fiberPolyData != nullptr) m_FiberPolyData = fiberPolyData; this->UpdateFiberGeometry(); this->GenerateFiberIds(); this->ColorFibersByOrientation(); } mitk::FiberBundle::~FiberBundle() { } mitk::FiberBundle::Pointer mitk::FiberBundle::GetDeepCopy() { mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(m_FiberPolyData); newFib->SetFiberColors(this->m_FiberColors); newFib->SetFiberWeights(this->m_FiberWeights); return newFib; } vtkSmartPointer mitk::FiberBundle::GeneratePolyDataByIds(std::vector fiberIds) { vtkSmartPointer newFiberPolyData = vtkSmartPointer::New(); vtkSmartPointer newLineSet = vtkSmartPointer::New(); vtkSmartPointer newPointSet = vtkSmartPointer::New(); auto finIt = fiberIds.begin(); while ( finIt != fiberIds.end() ) { if (*finIt < 0 || *finIt>GetNumFibers()){ MITK_INFO << "FiberID can not be negative or >NumFibers!!! check id Extraction!" << *finIt; break; } vtkSmartPointer fiber = m_FiberIdDataSet->GetCell(*finIt);//->DeepCopy(fiber); vtkSmartPointer fibPoints = fiber->GetPoints(); vtkSmartPointer newFiber = vtkSmartPointer::New(); newFiber->GetPointIds()->SetNumberOfIds( fibPoints->GetNumberOfPoints() ); for(int i=0; iGetNumberOfPoints(); i++) { newFiber->GetPointIds()->SetId(i, newPointSet->GetNumberOfPoints()); newPointSet->InsertNextPoint(fibPoints->GetPoint(i)[0], fibPoints->GetPoint(i)[1], fibPoints->GetPoint(i)[2]); } newLineSet->InsertNextCell(newFiber); ++finIt; } newFiberPolyData->SetPoints(newPointSet); newFiberPolyData->SetLines(newLineSet); return newFiberPolyData; } // merge two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::AddBundle(mitk::FiberBundle* fib) { if (fib==nullptr) { MITK_WARN << "trying to call AddBundle with nullptr argument"; return nullptr; } MITK_INFO << "Adding fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); // add current fiber bundle vtkSmartPointer weights = vtkSmartPointer::New(); weights->SetNumberOfValues(this->GetNumFibers()+fib->GetNumFibers()); unsigned int counter = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, this->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // add new fiber bundle for (int i=0; iGetFiberPolyData()->GetNumberOfCells(); i++) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p); vtkIdType id = vNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } weights->InsertValue(counter, fib->GetFiberWeight(i)); vNewLines->InsertNextCell(container); counter++; } // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(vNewPolyData); newFib->SetFiberWeights(weights); return newFib; } // subtract two fiber bundles mitk::FiberBundle::Pointer mitk::FiberBundle::SubtractBundle(mitk::FiberBundle* fib) { MITK_INFO << "Subtracting fibers"; vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = vtkSmartPointer::New(); vtkSmartPointer vNewPoints = vtkSmartPointer::New(); std::vector< std::vector< itk::Point > > points1; for( int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = GetItkPoint(points->GetPoint(0)); itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); points1.push_back( {start, end} ); } std::vector< std::vector< itk::Point > > points2; for( int i=0; iGetNumFibers(); i++ ) { vtkCell* cell = fib->GetFiberPolyData()->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; itk::Point start = GetItkPoint(points->GetPoint(0)); itk::Point end = GetItkPoint(points->GetPoint(numPoints-1)); points2.push_back( {start, end} ); } int progress = 0; std::vector< int > ids; #pragma omp parallel for for (int i=0; i<(int)points1.size(); i++) { #pragma omp critical { progress++; std::cout << (int)(100*(float)progress/points1.size()) << "%" << '\r'; cout.flush(); } bool match = false; for (unsigned int j=0; jGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (points==nullptr || numPoints<=0) continue; vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(points->GetPoint(j)); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } if(vNewLines->GetNumberOfCells()==0) return nullptr; // initialize PolyData vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); // initialize fiber bundle return mitk::FiberBundle::New(vNewPolyData); } itk::Point mitk::FiberBundle::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } /* * set PolyData (additional flag to recompute fiber geometry, default = true) */ void mitk::FiberBundle::SetFiberPolyData(vtkSmartPointer fiberPD, bool updateGeometry) { if (fiberPD == nullptr) this->m_FiberPolyData = vtkSmartPointer::New(); else m_FiberPolyData->DeepCopy(fiberPD); m_NumFibers = m_FiberPolyData->GetNumberOfLines(); if (updateGeometry) UpdateFiberGeometry(); GenerateFiberIds(); ColorFibersByOrientation(); } /* * return vtkPolyData */ vtkSmartPointer mitk::FiberBundle::GetFiberPolyData() const { return m_FiberPolyData; } void mitk::FiberBundle::ColorFibersByOrientation() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also PolyData needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= vtkPoints* extrPoints = nullptr; extrPoints = m_FiberPolyData->GetPoints(); int numOfPoints = 0; if (extrPoints!=nullptr) numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = 4; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(numOfPoints * componentSize); m_FiberColors->SetNumberOfComponents(componentSize); m_FiberColors->SetName("FIBER_COLORS"); int numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) return; /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); /* single fiber checkpoints: is number of points valid */ if (pointsPerFiber > 1) { /* operate on points of single fiber */ for (int i=0; i 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff1[2])); rgba[3] = (unsigned char) (255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = (unsigned char) (255.0 * std::fabs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::fabs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::fabs(diff2[2])); rgba[3] = (unsigned char) (255.0); } m_FiberColors->InsertTupleValue(idList[i], rgba); } } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; } else { MITK_DEBUG << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByCurvature(bool opacity, bool normalize) { double window = 5; //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = 4; m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * componentSize); m_FiberColors->SetNumberOfComponents(componentSize); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); vector< double > values; double min = 1; double max = 0; MITK_INFO << "Coloring fibers by curvature"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures for (int j=0; j > vectors; vnl_vector_fixed< float, 3 > meanV; meanV.fill(0.0); while(dist1) { double p1[3]; points->GetPoint(c-1, p1); double p2[3]; points->GetPoint(c, p2); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==j) meanV += v; c--; } c = j; dist = 0; while(distGetPoint(c, p1); double p2[3]; points->GetPoint(c+1, p2); vnl_vector_fixed< float, 3 > v; v[0] = p2[0]-p1[0]; v[1] = p2[1]-p1[1]; v[2] = p2[2]-p1[2]; dist += v.magnitude(); v.normalize(); vectors.push_back(v); if (c==j) meanV += v; c++; } meanV.normalize(); double dev = 0; for (unsigned int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); dev = 1.0-dev/180.0; values.push_back(dev); if (devmax) max = dev; } } unsigned int count = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); for (int j=0; j1) dev = 1; lookupTable->GetColor(dev, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); if (opacity) rgba[3] = (unsigned char) (255.0 * dev); else rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTupleValue(cell->GetPointId(j), rgba); count++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberOpacity(vtkDoubleArray* FAValArray) { for(long i=0; iGetNumberOfTuples(); i++) { double faValue = FAValArray->GetValue(i); faValue = faValue * 255.0; m_FiberColors->SetComponent(i,3, (unsigned char) faValue ); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ResetFiberOpacity() { for(long i=0; iGetNumberOfTuples(); i++) m_FiberColors->SetComponent(i,3, 255.0 ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByScalarMap(mitk::Image::Pointer FAimage, bool opacity, bool normalize) { mitkPixelTypeMultiplex3( ColorFibersByScalarMap, FAimage->GetPixelType(), FAimage, opacity, normalize ); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } template void mitk::FiberBundle::ColorFibersByScalarMap(const mitk::PixelType, mitk::Image::Pointer image, bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::ImagePixelReadAccessor readimage(image, image->GetVolumeData(0)); unsigned char rgba[4] = {0,0,0,0}; vtkPoints* pointSet = m_FiberPolyData->GetPoints(); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); double min = 9999999; double max = -9999999; for(long i=0; iGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; double pixelValue = readimage.GetPixelByWorldCoordinates(px); if (pixelValue>max) max = pixelValue; if (pixelValueGetNumberOfPoints(); ++i) { Point3D px; px[0] = pointSet->GetPoint(i)[0]; px[1] = pointSet->GetPoint(i)[1]; px[2] = pointSet->GetPoint(i)[2]; double pixelValue = readimage.GetPixelByWorldCoordinates(px); if (normalize) pixelValue = (pixelValue-min)/(max-min); else if (pixelValue>1) pixelValue = 1; double color[3]; lookupTable->GetColor(1-pixelValue, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); if (opacity) rgba[3] = (unsigned char) (255.0 * pixelValue); else rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTupleValue(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::ColorFibersByFiberWeights(bool opacity, bool normalize) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); mitk::LookupTable::Pointer mitkLookup = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(0.0, 0.8); lookupTable->Build(); mitkLookup->SetVtkLookupTable(lookupTable); mitkLookup->SetType(mitk::LookupTable::JET); unsigned char rgba[4] = {0,0,0,0}; unsigned int counter = 0; float max = -999999; float min = 999999; for (int i=0; iGetFiberWeight(i); if (weight>max) max = weight; if (weightGetCell(i); int numPoints = cell->GetNumberOfPoints(); double weight = this->GetFiberWeight(i); for (int j=0; j1) v = 1; double color[3]; lookupTable->GetColor(1-v, color); rgba[0] = (unsigned char) (255.0 * color[0]); rgba[1] = (unsigned char) (255.0 * color[1]); rgba[2] = (unsigned char) (255.0 * color[2]); if (opacity) rgba[3] = (unsigned char) (255.0 * v); else rgba[3] = (unsigned char) (255.0); m_FiberColors->InsertTupleValue(counter, rgba); counter++; } } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::SetFiberColors(float r, float g, float b, float alpha) { m_FiberColors = vtkSmartPointer::New(); m_FiberColors->Allocate(m_FiberPolyData->GetNumberOfPoints() * 4); m_FiberColors->SetNumberOfComponents(4); m_FiberColors->SetName("FIBER_COLORS"); unsigned char rgba[4] = {0,0,0,0}; for(long i=0; iGetNumberOfPoints(); ++i) { rgba[0] = (unsigned char) r; rgba[1] = (unsigned char) g; rgba[2] = (unsigned char) b; rgba[3] = (unsigned char) alpha; m_FiberColors->InsertTupleValue(i, rgba); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } void mitk::FiberBundle::GenerateFiberIds() { if (m_FiberPolyData == nullptr) return; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInputData(m_FiberPolyData); idFiberFilter->CellIdsOn(); // idFiberFilter->PointIdsOn(); // point id's are not needed idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); idFiberFilter->FieldDataOn(); idFiberFilter->Update(); m_FiberIdDataSet = idFiberFilter->GetOutput(); } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(ItkUcharImgType* mask, bool anyPoint, bool invert, bool bothEnds, float fraction) { vtkSmartPointer PolyData = m_FiberPolyData; if (anyPoint) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); fibCopy->ResampleLinear(minSpacing/5); PolyData = fibCopy->GetFiberPolyData(); } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Extracting fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cellOriginal = m_FiberPolyData->GetCell(i); int numPointsOriginal = cellOriginal->GetNumberOfPoints(); vtkPoints* pointsOriginal = cellOriginal->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); if (numPoints>1 && numPointsOriginal) { if (anyPoint) { int inside = 0; int outside = 0; if (!invert) { for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetLargestPossibleRegion().IsInside(idx) && mask->GetPixel(idx) != 0 ) { inside++; if (fraction==0) break; } else outside++; } float current_fraction = 0.0; if (inside+outside>0) current_fraction = (float)inside/(inside+outside); if (current_fraction>fraction) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else { bool includeFiber = true; for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) ) { inside++; includeFiber = false; break; } else outside++; } if (includeFiber) { for (int k=0; kGetPoint(k); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } else { double* start = pointsOriginal->GetPoint(0); itk::Point itkStart; itkStart[0] = start[0]; itkStart[1] = start[1]; itkStart[2] = start[2]; itk::Index<3> idxStart; mask->TransformPhysicalPointToIndex(itkStart, idxStart); double* end = pointsOriginal->GetPoint(numPointsOriginal-1); itk::Point itkEnd; itkEnd[0] = end[0]; itkEnd[1] = end[1]; itkEnd[2] = end[2]; itk::Index<3> idxEnd; mask->TransformPhysicalPointToIndex(itkEnd, idxEnd); if (invert) { if (bothEnds) { if ( mask->GetPixel(idxStart) == 0 && mask->GetPixel(idxEnd) == 0 ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else if ( mask->GetPixel(idxStart) == 0 || mask->GetPixel(idxEnd) == 0 ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else { if (bothEnds) { if ( mask->GetPixel(idxStart) != 0 && mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart) && mask->GetLargestPossibleRegion().IsInside(idxEnd) ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } else if ( (mask->GetPixel(idxStart) != 0 && mask->GetLargestPossibleRegion().IsInside(idxStart)) || (mask->GetPixel(idxEnd) != 0 && mask->GetLargestPossibleRegion().IsInside(idxEnd)) ) { for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } } } } } vtkNewCells->InsertNextCell(container); } if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); return mitk::FiberBundle::New(newPolyData); } mitk::FiberBundle::Pointer mitk::FiberBundle::RemoveFibersOutside(ItkUcharImgType* mask, bool invert) { float minSpacing = 1; if(mask->GetSpacing()[0]GetSpacing()[1] && mask->GetSpacing()[0]GetSpacing()[2]) minSpacing = mask->GetSpacing()[0]; else if (mask->GetSpacing()[1] < mask->GetSpacing()[2]) minSpacing = mask->GetSpacing()[1]; else minSpacing = mask->GetSpacing()[2]; mitk::FiberBundle::Pointer fibCopy = this->GetDeepCopy(); fibCopy->ResampleLinear(minSpacing/10); vtkSmartPointer PolyData =fibCopy->GetFiberPolyData(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Cutting fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); if (numPoints>1) { int newNumPoints = 0; for (int j=0; jGetPoint(j); itk::Point itkP; itkP[0] = p[0]; itkP[1] = p[1]; itkP[2] = p[2]; itk::Index<3> idx; mask->TransformPhysicalPointToIndex(itkP, idx); if ( mask->GetPixel(idx) != 0 && mask->GetLargestPossibleRegion().IsInside(idx) && !invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if ( (mask->GetPixel(idx) == 0 || !mask->GetLargestPossibleRegion().IsInside(idx)) && invert ) { vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); newNumPoints++; } else if (newNumPoints>0) { vtkNewCells->InsertNextCell(container); newNumPoints = 0; container = vtkSmartPointer::New(); } } if (newNumPoints>0) vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return nullptr; vtkSmartPointer newPolyData = vtkSmartPointer::New(); newPolyData->SetPoints(vtkNewPoints); newPolyData->SetLines(vtkNewCells); mitk::FiberBundle::Pointer newFib = mitk::FiberBundle::New(newPolyData); newFib->Compress(0.1); return newFib; } mitk::FiberBundle::Pointer mitk::FiberBundle::ExtractFiberSubset(DataNode* roi, DataStorage* storage) { if (roi==nullptr || !(dynamic_cast(roi->GetData()) || dynamic_cast(roi->GetData())) ) return nullptr; std::vector tmp = ExtractFiberIdSubset(roi, storage); if (tmp.size()<=0) return mitk::FiberBundle::New(); vtkSmartPointer pTmp = GeneratePolyDataByIds(tmp); return mitk::FiberBundle::New(pTmp); } std::vector mitk::FiberBundle::ExtractFiberIdSubset(DataNode *roi, DataStorage* storage) { std::vector result; if (roi==nullptr || roi->GetData()==nullptr) return result; mitk::PlanarFigureComposite::Pointer pfc = dynamic_cast(roi->GetData()); if (!pfc.IsNull()) // handle composite { DataStorage::SetOfObjects::ConstPointer children = storage->GetDerivations(roi); if (children->size()==0) return result; switch (pfc->getOperationType()) { case 0: // AND { MITK_INFO << "AND"; result = this->ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { std::vector inRoi = this->ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(std::min(result.size(),inRoi.size())); it = std::set_intersection(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( it - rest.begin() ); result = rest; } break; } case 1: // OR { MITK_INFO << "OR"; result = ExtractFiberIdSubset(children->ElementAt(0), storage); std::vector::iterator it; for (unsigned int i=1; iSize(); ++i) { it = result.end(); std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); result.insert(it, inRoi.begin(), inRoi.end()); } // remove duplicates sort(result.begin(), result.end()); it = unique(result.begin(), result.end()); result.resize( it - result.begin() ); break; } case 2: // NOT { MITK_INFO << "NOT"; for(long i=0; iGetNumFibers(); i++) result.push_back(i); std::vector::iterator it; for (unsigned int i=0; iSize(); ++i) { std::vector inRoi = ExtractFiberIdSubset(children->ElementAt(i), storage); std::vector rest(result.size()-inRoi.size()); it = std::set_difference(result.begin(), result.end(), inRoi.begin(), inRoi.end(), rest.begin() ); rest.resize( it - rest.begin() ); result = rest; } break; } } } else if ( dynamic_cast(roi->GetData()) ) // actual extraction { if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarPoly = dynamic_cast(roi->GetData()); //create vtkPolygon using controlpoints from planarFigure polygon vtkSmartPointer polygonVtk = vtkSmartPointer::New(); for (unsigned int i=0; iGetNumberOfControlPoints(); ++i) { itk::Point p = planarPoly->GetWorldControlPoint(i); vtkIdType id = polygonVtk->GetPoints()->InsertNextPoint(p[0], p[1], p[2] ); polygonVtk->GetPointIds()->InsertNextId(id); } MITK_INFO << "Extracting with polygon"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); double tolerance = 0.001; // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection double pcoords[3] = {0,0,0}; int subId = 0; int iD = polygonVtk->IntersectWithLine(p1, p2, tolerance, t, x, pcoords, subId); if (iD!=0) { result.push_back(i); break; } } } } else if ( dynamic_cast(roi->GetData()) ) { mitk::PlanarFigure::Pointer planarFigure = dynamic_cast(roi->GetData()); Vector3D planeNormal = planarFigure->GetPlaneGeometry()->GetNormal(); planeNormal.Normalize(); //calculate circle radius mitk::Point3D V1w = planarFigure->GetWorldControlPoint(0); //centerPoint mitk::Point3D V2w = planarFigure->GetWorldControlPoint(1); //radiusPoint double radius = V1w.EuclideanDistanceTo(V2w); radius *= radius; MITK_INFO << "Extracting with circle"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p1); double p2[3] = {0,0,0}; points->GetPoint(j+1, p2); // Outputs double t = 0; // Parametric coordinate of intersection (0 (corresponding to p1) to 1 (corresponding to p2)) double x[3] = {0,0,0}; // The coordinate of the intersection int iD = vtkPlane::IntersectWithLine(p1,p2,planeNormal.GetDataPointer(),V1w.GetDataPointer(),t,x); if (iD!=0) { double dist = (x[0]-V1w[0])*(x[0]-V1w[0])+(x[1]-V1w[1])*(x[1]-V1w[1])+(x[2]-V1w[2])*(x[2]-V1w[2]); if( dist <= radius) { result.push_back(i); break; } } } } } return result; } return result; } void mitk::FiberBundle::UpdateFiberGeometry() { vtkSmartPointer cleaner = vtkSmartPointer::New(); cleaner->SetInputData(m_FiberPolyData); cleaner->PointMergingOff(); cleaner->Update(); m_FiberPolyData = cleaner->GetOutput(); m_FiberLengths.clear(); m_MeanFiberLength = 0; m_MedianFiberLength = 0; m_LengthStDev = 0; m_NumFibers = m_FiberPolyData->GetNumberOfCells(); if (m_FiberColors==nullptr || m_FiberColors->GetNumberOfTuples()!=m_FiberPolyData->GetNumberOfPoints()) this->ColorFibersByOrientation(); if (m_FiberWeights->GetSize()!=m_NumFibers) { m_FiberWeights = vtkSmartPointer::New(); m_FiberWeights->SetName("FIBER_WEIGHTS"); m_FiberWeights->SetNumberOfValues(m_NumFibers); this->SetFiberWeights(1); } if (m_NumFibers<=0) // no fibers present; apply default geometry { m_MinFiberLength = 0; m_MaxFiberLength = 0; mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(false); float b[] = {0, 1, 0, 1, 0, 1}; geometry->SetFloatBounds(b); SetGeometry(geometry); return; } double b[6]; m_FiberPolyData->GetBounds(b); // calculate statistics for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); float length = 0; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); float dist = std::sqrt((p1[0]-p2[0])*(p1[0]-p2[0])+(p1[1]-p2[1])*(p1[1]-p2[1])+(p1[2]-p2[2])*(p1[2]-p2[2])); length += dist; } m_FiberLengths.push_back(length); m_MeanFiberLength += length; if (i==0) { m_MinFiberLength = length; m_MaxFiberLength = length; } else { if (lengthm_MaxFiberLength) m_MaxFiberLength = length; } } m_MeanFiberLength /= m_NumFibers; std::vector< float > sortedLengths = m_FiberLengths; std::sort(sortedLengths.begin(), sortedLengths.end()); for (int i=0; i1) m_LengthStDev /= (m_NumFibers-1); else m_LengthStDev = 0; m_LengthStDev = std::sqrt(m_LengthStDev); m_MedianFiberLength = sortedLengths.at(m_NumFibers/2); mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetFloatBounds(b); this->SetGeometry(geometry); m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } float mitk::FiberBundle::GetFiberWeight(unsigned int fiber) const { return m_FiberWeights->GetValue(fiber); } void mitk::FiberBundle::SetFiberWeights(float newWeight) { for (int i=0; iGetSize(); i++) m_FiberWeights->SetValue(i, newWeight); } void mitk::FiberBundle::SetFiberWeights(vtkSmartPointer weights) { if (m_NumFibers!=weights->GetSize()) { MITK_INFO << "Weights array not equal to number of fibers!"; return; } for (int i=0; iGetSize(); i++) m_FiberWeights->SetValue(i, weights->GetValue(i)); m_FiberWeights->SetName("FIBER_WEIGHTS"); } void mitk::FiberBundle::SetFiberWeight(unsigned int fiber, float weight) { m_FiberWeights->SetValue(fiber, weight); } void mitk::FiberBundle::SetFiberColors(vtkSmartPointer fiberColors) { for(long i=0; iGetNumberOfPoints(); ++i) { unsigned char source[4] = {0,0,0,0}; fiberColors->GetTupleValue(i, source); unsigned char target[4] = {0,0,0,0}; target[0] = source[0]; target[1] = source[1]; target[2] = source[2]; target[3] = source[3]; m_FiberColors->InsertTupleValue(i, target); } m_UpdateTime3D.Modified(); m_UpdateTime2D.Modified(); } itk::Matrix< double, 3, 3 > mitk::FiberBundle::TransformMatrix(itk::Matrix< double, 3, 3 > m, double rx, double ry, double rz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; itk::Matrix< double, 3, 3 > rotX; rotX.SetIdentity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; itk::Matrix< double, 3, 3 > rotY; rotY.SetIdentity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; itk::Matrix< double, 3, 3 > rotZ; rotZ.SetIdentity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; itk::Matrix< double, 3, 3 > rot = rotZ*rotY*rotX; m = rot*m; return m; } itk::Point mitk::FiberBundle::TransformPoint(vnl_vector_fixed< double, 3 > point, double rx, double ry, double rz, double tx, double ty, double tz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); point[0] -= center[0]; point[1] -= center[1]; point[2] -= center[2]; point = rot*point; point[0] += center[0]+tx; point[1] += center[1]+ty; point[2] += center[2]+tz; itk::Point out; out[0] = point[0]; out[1] = point[1]; out[2] = point[2]; return out; } void mitk::FiberBundle::TransformFibers(double rx, double ry, double rz, double tx, double ty, double tz) { rx = rx*M_PI/180; ry = ry*M_PI/180; rz = rz*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(rx); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(rx); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(ry); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(ry); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(rz); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(rz); rotZ[1][0] = -rotZ[0][1]; vnl_matrix_fixed< double, 3, 3 > rot = rotZ*rotY*rotX; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rot*dir; dir[0] += center[0]+tx; dir[1] += center[1]+ty; dir[2] += center[2]+tz; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RotateAroundAxis(double x, double y, double z) { x = x*M_PI/180; y = y*M_PI/180; z = z*M_PI/180; vnl_matrix_fixed< double, 3, 3 > rotX; rotX.set_identity(); rotX[1][1] = cos(x); rotX[2][2] = rotX[1][1]; rotX[1][2] = -sin(x); rotX[2][1] = -rotX[1][2]; vnl_matrix_fixed< double, 3, 3 > rotY; rotY.set_identity(); rotY[0][0] = cos(y); rotY[2][2] = rotY[0][0]; rotY[0][2] = sin(y); rotY[2][0] = -rotY[0][2]; vnl_matrix_fixed< double, 3, 3 > rotZ; rotZ.set_identity(); rotZ[0][0] = cos(z); rotZ[1][1] = rotZ[0][0]; rotZ[0][1] = -sin(z); rotZ[1][0] = -rotZ[0][1]; mitk::BaseGeometry::Pointer geom = this->GetGeometry(); mitk::Point3D center = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vnl_vector_fixed< double, 3 > dir; dir[0] = p[0]-center[0]; dir[1] = p[1]-center[1]; dir[2] = p[2]-center[2]; dir = rotZ*rotY*rotX*dir; dir[0] += center[0]; dir[1] += center[1]; dir[2] += center[2]; vtkIdType id = vtkNewPoints->InsertNextPoint(dir.data_block()); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ScaleFibers(double x, double y, double z, bool subtractCenter) { MITK_INFO << "Scaling fibers"; boost::progress_display disp(m_NumFibers); mitk::BaseGeometry* geom = this->GetGeometry(); mitk::Point3D c = geom->GetCenter(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); if (subtractCenter) { p[0] -= c[0]; p[1] -= c[1]; p[2] -= c[2]; } p[0] *= x; p[1] *= y; p[2] *= z; if (subtractCenter) { p[0] += c[0]; p[1] += c[1]; p[2] += c[2]; } vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::TranslateFibers(double x, double y, double z) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[0] += x; p[1] += y; p[2] += z; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::MirrorFibers(unsigned int axis) { if (axis>2) return; MITK_INFO << "Mirroring fibers"; boost::progress_display disp(m_NumFibers); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); p[axis] = -p[axis]; vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::RemoveDir(vnl_vector_fixed dir, double threshold) { dir.normalize(); vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); bool discard = false; for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); vnl_vector_fixed< double, 3 > v1; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; if (v1.magnitude()>0.001) { v1.normalize(); if (fabs(dot_product(v1,dir))>threshold) { discard = true; break; } } } if (!discard) { for (int j=0; jGetPoint(j, p1); vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); // UpdateColorCoding(); // UpdateFiberGeometry(); } bool mitk::FiberBundle::ApplyCurvatureThreshold(float minRadius, bool deleteFibers) { if (minRadius<0) return true; vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Applying curvature threshold"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); for (int i=0; iGetNumberOfCells(); i++) { ++disp ; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); // calculate curvatures vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j, p1); double p2[3]; points->GetPoint(j+1, p2); double p3[3]; points->GetPoint(j+2, p3); vnl_vector_fixed< float, 3 > v1, v2, v3; v1[0] = p2[0]-p1[0]; v1[1] = p2[1]-p1[1]; v1[2] = p2[2]-p1[2]; v2[0] = p3[0]-p2[0]; v2[1] = p3[1]-p2[1]; v2[2] = p3[2]-p2[2]; v3[0] = p1[0]-p3[0]; v3[1] = p1[1]-p3[1]; v3[2] = p1[2]-p3[2]; float a = v1.magnitude(); float b = v2.magnitude(); float c = v3.magnitude(); float r = a*b*c/std::sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c)); // radius of triangle via Heron's formula (area of triangle) vtkIdType id = vtkNewPoints->InsertNextPoint(p1); container->GetPointIds()->InsertNextId(id); if (deleteFibers && rInsertNextCell(container); container = vtkSmartPointer::New(); } else if (j==numPoints-3) { id = vtkNewPoints->InsertNextPoint(p2); container->GetPointIds()->InsertNextId(id); id = vtkNewPoints->InsertNextPoint(p3); container->GetPointIds()->InsertNextId(id); vtkNewCells->InsertNextCell(container); } } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveShortFibers(float lengthInMM) { MITK_INFO << "Removing short fibers"; if (lengthInMM<=0 || lengthInMMm_MaxFiberLength) // can't remove all fibers { MITK_WARN << "Process aborted. No fibers would be left!"; return false; } vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); float min = m_MaxFiberLength; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)>=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); if (m_FiberLengths.at(i)GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } bool mitk::FiberBundle::RemoveLongFibers(float lengthInMM) { if (lengthInMM<=0 || lengthInMM>m_MaxFiberLength) return true; if (lengthInMM vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Removing long fibers"; boost::progress_display disp(m_NumFibers); for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); if (m_FiberLengths.at(i)<=lengthInMM) { vtkSmartPointer container = vtkSmartPointer::New(); for (int j=0; jGetPoint(j); vtkIdType id = vtkNewPoints->InsertNextPoint(p); container->GetPointIds()->InsertNextId(id); } vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()<=0) return false; m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); return true; } void mitk::FiberBundle::ResampleSpline(float pointDistance, double tension, double continuity, double bias ) { if (pointDistance<=0) return; vtkSmartPointer vtkSmoothPoints = vtkSmartPointer::New(); //in smoothpoints the interpolated points representing a fiber are stored. //in vtkcells all polylines are stored, actually all id's of them are stored vtkSmartPointer vtkSmoothCells = vtkSmartPointer::New(); //cellcontainer for smoothed lines vtkIdType pointHelperCnt = 0; MITK_INFO << "Smoothing fibers"; vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); boost::progress_display disp(m_NumFibers); #pragma omp parallel for for (int i=0; i newPoints = vtkSmartPointer::New(); float length = 0; float weight = 1; #pragma omp critical { length = m_FiberLengths.at(i); weight = m_FiberWeights->GetValue(i); ++disp; vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jInsertNextPoint(points->GetPoint(j)); } int sampling = std::ceil(length/pointDistance); vtkSmartPointer xSpline = vtkSmartPointer::New(); vtkSmartPointer ySpline = vtkSmartPointer::New(); vtkSmartPointer zSpline = vtkSmartPointer::New(); xSpline->SetDefaultBias(bias); xSpline->SetDefaultTension(tension); xSpline->SetDefaultContinuity(continuity); ySpline->SetDefaultBias(bias); ySpline->SetDefaultTension(tension); ySpline->SetDefaultContinuity(continuity); zSpline->SetDefaultBias(bias); zSpline->SetDefaultTension(tension); zSpline->SetDefaultContinuity(continuity); vtkSmartPointer spline = vtkSmartPointer::New(); spline->SetXSpline(xSpline); spline->SetYSpline(ySpline); spline->SetZSpline(zSpline); spline->SetPoints(newPoints); vtkSmartPointer functionSource = vtkSmartPointer::New(); functionSource->SetParametricFunction(spline); functionSource->SetUResolution(sampling); functionSource->SetVResolution(sampling); functionSource->SetWResolution(sampling); functionSource->Update(); vtkPolyData* outputFunction = functionSource->GetOutput(); vtkPoints* tmpSmoothPnts = outputFunction->GetPoints(); //smoothPoints of current fiber vtkSmartPointer smoothLine = vtkSmartPointer::New(); smoothLine->GetPointIds()->SetNumberOfIds(tmpSmoothPnts->GetNumberOfPoints()); #pragma omp critical { for (int j=0; jGetNumberOfPoints(); j++) { smoothLine->GetPointIds()->SetId(j, j+pointHelperCnt); vtkSmoothPoints->InsertNextPoint(tmpSmoothPnts->GetPoint(j)); } newFiberWeights->SetValue(vtkSmoothCells->GetNumberOfCells(), weight); vtkSmoothCells->InsertNextCell(smoothLine); pointHelperCnt += tmpSmoothPnts->GetNumberOfPoints(); } } SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkSmoothPoints); m_FiberPolyData->SetLines(vtkSmoothCells); this->SetFiberPolyData(m_FiberPolyData, true); } void mitk::FiberBundle::ResampleSpline(float pointDistance) { ResampleSpline(pointDistance, 0, 0, 0 ); } unsigned long mitk::FiberBundle::GetNumberOfPoints() const { unsigned long points = 0; for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); points += cell->GetNumberOfPoints(); } return points; } void mitk::FiberBundle::Compress(float error) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Compressing fibers"; unsigned long numRemovedPoints = 0; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } // calculate curvatures int numPoints = vertices.size(); std::vector< int > removedPoints; removedPoints.resize(numPoints, 0); removedPoints[0]=-1; removedPoints[numPoints-1]=-1; vtkSmartPointer container = vtkSmartPointer::New(); int remCounter = 0; bool pointFound = true; while (pointFound) { pointFound = false; double minError = error; int removeIndex = -1; for (unsigned int j=0; j candV = vertices.at(j); int validP = -1; vnl_vector_fixed< double, 3 > pred; for (int k=j-1; k>=0; k--) if (removedPoints[k]<=0) { pred = vertices.at(k); validP = k; break; } int validS = -1; vnl_vector_fixed< double, 3 > succ; for (int k=j+1; k=0 && validS>=0) { double a = (candV-pred).magnitude(); double b = (candV-succ).magnitude(); double c = (pred-succ).magnitude(); double s=0.5*(a+b+c); double hc=(2.0/c)*sqrt(fabs(s*(s-a)*(s-b)*(s-c))); if (hcInsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); numRemovedPoints += remCounter; vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { MITK_INFO << "Removed points: " << numRemovedPoints; SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } void mitk::FiberBundle::ResampleLinear(double pointDistance) { vtkSmartPointer vtkNewPoints = vtkSmartPointer::New(); vtkSmartPointer vtkNewCells = vtkSmartPointer::New(); MITK_INFO << "Resampling fibers (linear)"; boost::progress_display disp(m_FiberPolyData->GetNumberOfCells()); vtkSmartPointer newFiberWeights = vtkSmartPointer::New(); newFiberWeights->SetName("FIBER_WEIGHTS"); newFiberWeights->SetNumberOfValues(m_NumFibers); #pragma omp parallel for for (int i=0; iGetNumberOfCells(); i++) { std::vector< vnl_vector_fixed< double, 3 > > vertices; float weight = 1; #pragma omp critical { ++disp; weight = m_FiberWeights->GetValue(i); vtkCell* cell = m_FiberPolyData->GetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, cand); vnl_vector_fixed< double, 3 > candV; candV[0]=cand[0]; candV[1]=cand[1]; candV[2]=cand[2]; vertices.push_back(candV); } } vtkSmartPointer container = vtkSmartPointer::New(); vnl_vector_fixed< double, 3 > lastV = vertices.at(0); #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(lastV.data_block()); container->GetPointIds()->InsertNextId(id); } for (unsigned int j=1; j vec = vertices.at(j) - lastV; double new_dist = vec.magnitude(); if (new_dist >= pointDistance) { vnl_vector_fixed< double, 3 > newV = lastV; if ( new_dist-pointDistance <= mitk::eps ) { vec.normalize(); newV += vec * pointDistance; } else { // intersection between sphere (radius 'pointDistance', center 'lastV') and line (direction 'd' and point 'p') vnl_vector_fixed< double, 3 > p = vertices.at(j-1); vnl_vector_fixed< double, 3 > d = vertices.at(j) - p; double a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]; double b = 2 * (d[0] * (p[0] - lastV[0]) + d[1] * (p[1] - lastV[1]) + d[2] * (p[2] - lastV[2])); double c = (p[0] - lastV[0])*(p[0] - lastV[0]) + (p[1] - lastV[1])*(p[1] - lastV[1]) + (p[2] - lastV[2])*(p[2] - lastV[2]) - pointDistance*pointDistance; double v1 =(-b + std::sqrt(b*b-4*a*c))/(2*a); double v2 =(-b - std::sqrt(b*b-4*a*c))/(2*a); if (v1>0) newV = p + d * v1; else if (v2>0) newV = p + d * v2; else MITK_INFO << "ERROR1 - linear resampling"; j--; } #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(newV.data_block()); container->GetPointIds()->InsertNextId(id); } lastV = newV; } else if (j==vertices.size()-1 && new_dist>0.0001) { #pragma omp critical { vtkIdType id = vtkNewPoints->InsertNextPoint(vertices.at(j).data_block()); container->GetPointIds()->InsertNextId(id); } } } #pragma omp critical { newFiberWeights->SetValue(vtkNewCells->GetNumberOfCells(), weight); vtkNewCells->InsertNextCell(container); } } if (vtkNewCells->GetNumberOfCells()>0) { SetFiberWeights(newFiberWeights); m_FiberPolyData = vtkSmartPointer::New(); m_FiberPolyData->SetPoints(vtkNewPoints); m_FiberPolyData->SetLines(vtkNewCells); this->SetFiberPolyData(m_FiberPolyData, true); } } // reapply selected colorcoding in case PolyData structure has changed bool mitk::FiberBundle::Equals(mitk::FiberBundle* fib, double eps) { if (fib==nullptr) { MITK_INFO << "Reference bundle is nullptr!"; return false; } if (m_NumFibers!=fib->GetNumFibers()) { MITK_INFO << "Unequal number of fibers!"; MITK_INFO << m_NumFibers << " vs. " << fib->GetNumFibers(); return false; } for (int i=0; iGetCell(i); int numPoints = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); vtkCell* cell2 = fib->GetFiberPolyData()->GetCell(i); int numPoints2 = cell2->GetNumberOfPoints(); vtkPoints* points2 = cell2->GetPoints(); if (numPoints2!=numPoints) { MITK_INFO << "Unequal number of points in fiber " << i << "!"; MITK_INFO << numPoints2 << " vs. " << numPoints; return false; } for (int j=0; jGetPoint(j); double* p2 = points2->GetPoint(j); if (fabs(p1[0]-p2[0])>eps || fabs(p1[1]-p2[1])>eps || fabs(p1[2]-p2[2])>eps) { MITK_INFO << "Unequal points in fiber " << i << " at position " << j << "!"; MITK_INFO << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2]; MITK_INFO << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2]; return false; } } } return true; } void mitk::FiberBundle::PrintSelf(std::ostream &os, itk::Indent indent) const { os << indent << this->GetNameOfClass() << ":\n"; os << indent << "Number of fibers: " << this->GetNumFibers() << std::endl; os << indent << "Min. fiber length: " << this->GetMinFiberLength() << std::endl; os << indent << "Max. fiber length: " << this->GetMaxFiberLength() << std::endl; os << indent << "Mean fiber length: " << this->GetMeanFiberLength() << std::endl; os << indent << "Median fiber length: " << this->GetMedianFiberLength() << std::endl; os << indent << "STDEV fiber length: " << this->GetLengthStDev() << std::endl; os << indent << "Number of points: " << this->GetNumberOfPoints() << std::endl; - os << indent << "" << this->GetNumberOfPoints() << std::endl; Superclass::PrintSelf(os, indent); } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundle::UpdateOutputInformation() { } void mitk::FiberBundle::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundle::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundle::VerifyRequestedRegion() { return true; } void mitk::FiberBundle::SetRequestedRegion(const itk::DataObject* ) { } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp index cf9a0b6af8..887f3f7a6b 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingView.cpp @@ -1,635 +1,635 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkGibbsTrackingView.h" // Qt #include #include #include #include // MITK #include #include #include #include #include // ITK #include #include #include // MISC #include QmitkTrackingWorker::QmitkTrackingWorker(QmitkGibbsTrackingView* view) - : m_View(view) + : m_View(view) { } void QmitkTrackingWorker::run() { - m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); - - m_View->m_GlobalTracker->SetQBallImage(m_View->m_ItkQBallImage); - m_View->m_GlobalTracker->SetTensorImage(m_View->m_ItkTensorImage); - m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); - m_View->m_GlobalTracker->SetStartTemperature((float)m_View->m_Controls->m_StartTempSlider->value()/100); - m_View->m_GlobalTracker->SetEndTemperature((float)m_View->m_Controls->m_EndTempSlider->value()/10000); - m_View->m_GlobalTracker->SetIterations(m_View->m_Controls->m_IterationsBox->text().toDouble()); - m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); - m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); - m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); - m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); - m_View->m_GlobalTracker->SetMinFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); - m_View->m_GlobalTracker->SetCurvatureThreshold(cos((float)m_View->m_Controls->m_CurvatureThresholdSlider->value()*M_PI/180)); - m_View->m_GlobalTracker->SetRandomSeed(m_View->m_Controls->m_RandomSeedSlider->value()); - try{ - m_View->m_GlobalTracker->Update(); - } - catch( mitk::Exception e ) - { - MITK_ERROR << "Internal error occured: " << e.what() << "\nAborting"; - } - m_View->m_TrackingThread.quit(); + m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); + + m_View->m_GlobalTracker->SetQBallImage(m_View->m_ItkQBallImage); + m_View->m_GlobalTracker->SetTensorImage(m_View->m_ItkTensorImage); + m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); + m_View->m_GlobalTracker->SetStartTemperature((float)m_View->m_Controls->m_StartTempSlider->value()/100); + m_View->m_GlobalTracker->SetEndTemperature((float)m_View->m_Controls->m_EndTempSlider->value()/10000); + m_View->m_GlobalTracker->SetIterations(m_View->m_Controls->m_IterationsBox->text().toDouble()); + m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); + m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); + m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); + m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); + m_View->m_GlobalTracker->SetMinFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); + m_View->m_GlobalTracker->SetCurvatureThreshold(cos((float)m_View->m_Controls->m_CurvatureThresholdSlider->value()*M_PI/180)); + m_View->m_GlobalTracker->SetRandomSeed(m_View->m_Controls->m_RandomSeedSlider->value()); + try{ + m_View->m_GlobalTracker->Update(); + } + catch( mitk::Exception e ) + { + MITK_ERROR << "Internal error occured: " << e.what() << "\nAborting"; + } + m_View->m_TrackingThread.quit(); } const std::string QmitkGibbsTrackingView::VIEW_ID = "org.mitk.views.gibbstracking"; QmitkGibbsTrackingView::QmitkGibbsTrackingView() - : QmitkAbstractView() - , m_Controls( 0 ) - , m_TrackingNode(nullptr) - , m_FiberBundle(nullptr) - , m_MaskImage(nullptr) - , m_TensorImage(nullptr) - , m_QBallImage(nullptr) - , m_ItkQBallImage(nullptr) - , m_ItkTensorImage(nullptr) - , m_ImageNode(nullptr) - , m_MaskImageNode(nullptr) - , m_FiberBundleNode(nullptr) - , m_ThreadIsRunning(false) - , m_ElapsedTime(0) - , m_GlobalTracker(nullptr) - , m_TrackingWorker(this) + : QmitkAbstractView() + , m_Controls( 0 ) + , m_TrackingNode(nullptr) + , m_FiberBundle(nullptr) + , m_MaskImage(nullptr) + , m_TensorImage(nullptr) + , m_QBallImage(nullptr) + , m_ItkQBallImage(nullptr) + , m_ItkTensorImage(nullptr) + , m_ImageNode(nullptr) + , m_MaskImageNode(nullptr) + , m_FiberBundleNode(nullptr) + , m_ThreadIsRunning(false) + , m_ElapsedTime(0) + , m_GlobalTracker(nullptr) + , m_TrackingWorker(this) { - m_TrackingWorker.moveToThread(&m_TrackingThread); - connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); - connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); - connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); - m_TrackingTimer = new QTimer(this); + m_TrackingWorker.moveToThread(&m_TrackingThread); + connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); + connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); + connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); + m_TrackingTimer = new QTimer(this); } QmitkGibbsTrackingView::~QmitkGibbsTrackingView() { if (m_GlobalTracker.IsNull()) - return; + return; m_GlobalTracker->SetAbortTracking(true); m_TrackingThread.wait(); } // update tracking status and generate fiber bundle void QmitkGibbsTrackingView::TimerUpdate() { - UpdateTrackingStatus(); - GenerateFiberBundle(); + UpdateTrackingStatus(); + GenerateFiberBundle(); } // tell global tractography filter to stop after current step void QmitkGibbsTrackingView::StopGibbsTracking() { - if (m_GlobalTracker.IsNull()) - return; + if (m_GlobalTracker.IsNull()) + return; - m_GlobalTracker->SetAbortTracking(true); - m_Controls->m_TrackingStop->setEnabled(false); - m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); - m_TrackingNode = nullptr; + m_GlobalTracker->SetAbortTracking(true); + m_Controls->m_TrackingStop->setEnabled(false); + m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); + m_TrackingNode = nullptr; } // update gui elements and generate fiber bundle after tracking is finished void QmitkGibbsTrackingView::AfterThread() { - m_ThreadIsRunning = false; - m_TrackingTimer->stop(); + m_ThreadIsRunning = false; + m_TrackingTimer->stop(); - UpdateGUI(); + UpdateGUI(); - if( !m_GlobalTracker->GetIsInValidState() ) - { - QMessageBox::critical( nullptr, "Gibbs Tracking", "An internal error occured. Tracking aborted.\n Please check the log for details." ); - m_FiberBundleNode = nullptr; - return; - } - UpdateTrackingStatus(); + if( !m_GlobalTracker->GetIsInValidState() ) + { + QMessageBox::critical( nullptr, "Gibbs Tracking", "An internal error occured. Tracking aborted.\n Please check the log for details." ); + m_FiberBundleNode = nullptr; + return; + } + UpdateTrackingStatus(); - if(m_Controls->m_ParticleWeightSlider->value()==0) - { - m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); - m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); - } - if(m_Controls->m_ParticleWidthSlider->value()==0) - { - m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); - m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); - } - if(m_Controls->m_ParticleLengthSlider->value()==0) - { - m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); - m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); - } + if(m_Controls->m_ParticleWeightSlider->value()==0) + { + m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); + m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); + } + if(m_Controls->m_ParticleWidthSlider->value()==0) + { + m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); + m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); + } + if(m_Controls->m_ParticleLengthSlider->value()==0) + { + m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); + m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); + } - GenerateFiberBundle(); - m_FiberBundleNode = 0; - m_GlobalTracker = 0; + GenerateFiberBundle(); + m_FiberBundleNode = 0; + m_GlobalTracker = 0; - // images not needed anymore ( relevant only for computation ) - // we need to release them to remove the memory access block created through CastToItk<> calls - this->m_ItkQBallImage = 0; - this->m_ItkTensorImage = 0; + // images not needed anymore ( relevant only for computation ) + // we need to release them to remove the memory access block created through CastToItk<> calls + this->m_ItkQBallImage = 0; + this->m_ItkTensorImage = 0; } // start tracking timer and update gui elements before tracking is started void QmitkGibbsTrackingView::BeforeThread() { - m_ThreadIsRunning = true; - m_TrackingTime = QTime::currentTime(); - m_ElapsedTime = 0; - m_TrackingTimer->start(1000); + m_ThreadIsRunning = true; + m_TrackingTime = QTime::currentTime(); + m_ElapsedTime = 0; + m_TrackingTimer->start(1000); - UpdateGUI(); + UpdateGUI(); } // setup gui elements and signal/slot connections void QmitkGibbsTrackingView::CreateQtPartControl( QWidget *parent ) { - // build up qt view, unless already done - if ( !m_Controls ) - { - // create GUI widgets from the Qt Designer's .ui file - m_Controls = new Ui::QmitkGibbsTrackingViewControls; - m_Controls->setupUi( parent ); - - AdvancedSettings(); - - connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); - connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); - connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); - connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); - connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); - connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); - connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); - connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); - connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); - connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); - connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); - connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); - connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); - connect( m_Controls->m_CurvatureThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(SetCurvatureThreshold(int)) ); - connect( m_Controls->m_RandomSeedSlider, SIGNAL(valueChanged(int)), this, SLOT(SetRandomSeed(int)) ); - connect( m_Controls->m_OutputFileButton, SIGNAL(clicked()), this, SLOT(SetOutputFile()) ); - } + // build up qt view, unless already done + if ( !m_Controls ) + { + // create GUI widgets from the Qt Designer's .ui file + m_Controls = new Ui::QmitkGibbsTrackingViewControls; + m_Controls->setupUi( parent ); + + AdvancedSettings(); + + connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); + connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); + connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); + connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); + connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); + connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); + connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); + connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); + connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); + connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); + connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); + connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); + connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); + connect( m_Controls->m_CurvatureThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(SetCurvatureThreshold(int)) ); + connect( m_Controls->m_RandomSeedSlider, SIGNAL(valueChanged(int)), this, SLOT(SetRandomSeed(int)) ); + connect( m_Controls->m_OutputFileButton, SIGNAL(clicked()), this, SLOT(SetOutputFile()) ); + } } void QmitkGibbsTrackingView::SetFocus() { m_Controls->m_TrackingStart->setFocus(); } void QmitkGibbsTrackingView::SetInExBalance(int value) { - m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); + m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); } void QmitkGibbsTrackingView::SetFiberLength(int value) { - m_Controls->m_FiberLengthLabel->setText(QString::number(value)+"mm"); + m_Controls->m_FiberLengthLabel->setText(QString::number(value)+"mm"); } void QmitkGibbsTrackingView::SetRandomSeed(int value) { - if (value>=0) - m_Controls->m_RandomSeedLabel->setText(QString::number(value)); - else - m_Controls->m_RandomSeedLabel->setText("auto"); + if (value>=0) + m_Controls->m_RandomSeedLabel->setText(QString::number(value)); + else + m_Controls->m_RandomSeedLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleWeight(int value) { - if (value>0) - m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); - else - m_Controls->m_ParticleWeightLabel->setText("auto"); + if (value>0) + m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); + else + m_Controls->m_ParticleWeightLabel->setText("auto"); } void QmitkGibbsTrackingView::SetStartTemp(int value) { - m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); + m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); } void QmitkGibbsTrackingView::SetEndTemp(int value) { - m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); + m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); } void QmitkGibbsTrackingView::SetParticleWidth(int value) { - if (value>0) - m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); - else - m_Controls->m_ParticleWidthLabel->setText("auto"); + if (value>0) + m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); + else + m_Controls->m_ParticleWidthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleLength(int value) { - if (value>0) - m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); - else - m_Controls->m_ParticleLengthLabel->setText("auto"); + if (value>0) + m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); + else + m_Controls->m_ParticleLengthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetCurvatureThreshold(int value) { - m_Controls->m_CurvatureThresholdLabel->setText(QString::number(value)+"°"); + m_Controls->m_CurvatureThresholdLabel->setText(QString::number(value)+"°"); } // called if datamanager selection changes void QmitkGibbsTrackingView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { - if (m_ThreadIsRunning) - return; + if (m_ThreadIsRunning) + return; - m_ImageNode = nullptr; - m_MaskImageNode = nullptr; + m_ImageNode = nullptr; + m_MaskImageNode = nullptr; - // iterate all selected objects - for (mitk::DataNode::Pointer node: nodes) + // iterate all selected objects + for (mitk::DataNode::Pointer node: nodes) + { + if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + m_ImageNode = node; + else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + m_ImageNode = node; + else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { - if( node.IsNotNull() && dynamic_cast(node->GetData()) ) - m_ImageNode = node; - else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) - m_ImageNode = node; - else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) - { - mitk::Image::Pointer img = dynamic_cast(node->GetData()); - if (img->GetPixelType().GetPixelType()==itk::ImageIOBase::SCALAR) - m_MaskImageNode = node; - } + mitk::Image::Pointer img = dynamic_cast(node->GetData()); + if (img->GetPixelType().GetPixelType()==itk::ImageIOBase::SCALAR) + m_MaskImageNode = node; } + } - UpdateGUI(); + UpdateGUI(); } void QmitkGibbsTrackingView::NodeRemoved(const mitk::DataNode * node) { if (m_ThreadIsRunning) { if (node==m_TrackingNode.GetPointer()) { StopGibbsTracking(); } } } // update gui elements displaying trackings status void QmitkGibbsTrackingView::UpdateTrackingStatus() { - if (m_GlobalTracker.IsNull()) - return; - - m_ElapsedTime += m_TrackingTime.elapsed()/1000; - m_TrackingTime.restart(); - unsigned long hours = m_ElapsedTime/3600; - unsigned long minutes = (m_ElapsedTime%3600)/60; - unsigned long seconds = m_ElapsedTime%60; - - m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); - - m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); - m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); - m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); - m_Controls->m_CurrentStepLabel->setText( QString::number(100*m_GlobalTracker->GetCurrentIteration()/m_GlobalTracker->GetIterations())+"%" ); - m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); + if (m_GlobalTracker.IsNull()) + return; + + m_ElapsedTime += m_TrackingTime.elapsed()/1000; + m_TrackingTime.restart(); + unsigned long hours = m_ElapsedTime/3600; + unsigned long minutes = (m_ElapsedTime%3600)/60; + unsigned long seconds = m_ElapsedTime%60; + + m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); + + m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); + m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); + m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); + m_Controls->m_CurrentStepLabel->setText( QString::number(100*m_GlobalTracker->GetCurrentIteration()/m_GlobalTracker->GetIterations())+"%" ); + m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); } // update gui elements (enable/disable elements and set tooltips) void QmitkGibbsTrackingView::UpdateGUI() { - if (m_ImageNode.IsNotNull()) - { - m_Controls->m_QballImageLabel->setText(m_ImageNode->GetName().c_str()); - m_Controls->m_DataFrame->setTitle("Input Data"); - } - else - { - m_Controls->m_QballImageLabel->setText("mandatory"); - m_Controls->m_DataFrame->setTitle("Please Select Input Data"); - } - if (m_MaskImageNode.IsNotNull()) - m_Controls->m_MaskImageLabel->setText(m_MaskImageNode->GetName().c_str()); - else - m_Controls->m_MaskImageLabel->setText("optional"); + if (m_ImageNode.IsNotNull()) + { + m_Controls->m_QballImageLabel->setText(m_ImageNode->GetName().c_str()); + m_Controls->m_DataFrame->setTitle("Input Data"); + } + else + { + m_Controls->m_QballImageLabel->setText("mandatory"); + m_Controls->m_DataFrame->setTitle("Please Select Input Data"); + } + if (m_MaskImageNode.IsNotNull()) + m_Controls->m_MaskImageLabel->setText(m_MaskImageNode->GetName().c_str()); + else + m_Controls->m_MaskImageLabel->setText("optional"); - if (!m_ThreadIsRunning && m_ImageNode.IsNotNull()) - { - m_Controls->m_TrackingStop->setEnabled(false); - m_Controls->m_TrackingStart->setEnabled(true); - m_Controls->m_LoadTrackingParameters->setEnabled(true); - m_Controls->m_IterationsBox->setEnabled(true); - m_Controls->m_AdvancedFrame->setEnabled(true); - m_Controls->m_TrackingStop->setText("Stop Tractography"); - m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); - m_Controls->m_TrackingStop->setToolTip(""); - } - else if (!m_ThreadIsRunning) - { - m_Controls->m_TrackingStop->setEnabled(false); - m_Controls->m_TrackingStart->setEnabled(false); - m_Controls->m_LoadTrackingParameters->setEnabled(true); - m_Controls->m_IterationsBox->setEnabled(true); - m_Controls->m_AdvancedFrame->setEnabled(true); - m_Controls->m_TrackingStop->setText("Stop Tractography"); - m_Controls->m_TrackingStart->setToolTip("No Q-Ball image selected."); - m_Controls->m_TrackingStop->setToolTip(""); - } - else - { - m_Controls->m_TrackingStop->setEnabled(true); - m_Controls->m_TrackingStart->setEnabled(false); - m_Controls->m_LoadTrackingParameters->setEnabled(false); - m_Controls->m_IterationsBox->setEnabled(false); - m_Controls->m_AdvancedFrame->setEnabled(false); - m_Controls->m_AdvancedFrame->setVisible(false); - m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); - m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); - m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); - } + if (!m_ThreadIsRunning && m_ImageNode.IsNotNull()) + { + m_Controls->m_TrackingStop->setEnabled(false); + m_Controls->m_TrackingStart->setEnabled(true); + m_Controls->m_LoadTrackingParameters->setEnabled(true); + m_Controls->m_IterationsBox->setEnabled(true); + m_Controls->m_AdvancedFrame->setEnabled(true); + m_Controls->m_TrackingStop->setText("Stop Tractography"); + m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); + m_Controls->m_TrackingStop->setToolTip(""); + } + else if (!m_ThreadIsRunning) + { + m_Controls->m_TrackingStop->setEnabled(false); + m_Controls->m_TrackingStart->setEnabled(false); + m_Controls->m_LoadTrackingParameters->setEnabled(true); + m_Controls->m_IterationsBox->setEnabled(true); + m_Controls->m_AdvancedFrame->setEnabled(true); + m_Controls->m_TrackingStop->setText("Stop Tractography"); + m_Controls->m_TrackingStart->setToolTip("No Q-Ball image selected."); + m_Controls->m_TrackingStop->setToolTip(""); + } + else + { + m_Controls->m_TrackingStop->setEnabled(true); + m_Controls->m_TrackingStart->setEnabled(false); + m_Controls->m_LoadTrackingParameters->setEnabled(false); + m_Controls->m_IterationsBox->setEnabled(false); + m_Controls->m_AdvancedFrame->setEnabled(false); + m_Controls->m_AdvancedFrame->setVisible(false); + m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); + m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); + m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); + } } // show/hide advanced settings frame void QmitkGibbsTrackingView::AdvancedSettings() { - m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); + m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); } // set mask image data node void QmitkGibbsTrackingView::SetMask() { - QList nodes = GetDataManagerSelection(); - if (nodes.empty()) - { - m_MaskImageNode = nullptr; - m_Controls->m_MaskImageLabel->setText("-"); - return; - } + QList nodes = GetDataManagerSelection(); + if (nodes.empty()) + { + m_MaskImageNode = nullptr; + m_Controls->m_MaskImageLabel->setText("-"); + return; + } - for (auto node: nodes) + for (auto node: nodes) + { + if (node.IsNotNull() && dynamic_cast(node->GetData())) { - if (node.IsNotNull() && dynamic_cast(node->GetData())) - { - m_MaskImageNode = node; - m_Controls->m_MaskImageLabel->setText(node->GetName().c_str()); - return; - } + m_MaskImageNode = node; + m_Controls->m_MaskImageLabel->setText(node->GetName().c_str()); + return; } + } } // check for mask and qbi and start tracking thread void QmitkGibbsTrackingView::StartGibbsTracking() { - if(m_ThreadIsRunning) - { - MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; - return; - } - m_GlobalTracker = nullptr; + if(m_ThreadIsRunning) + { + MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; + return; + } + m_GlobalTracker = nullptr; - if (m_ImageNode.IsNull()) - { - QMessageBox::information( nullptr, "Warning", "Please load and select a qball image before starting image processing."); - return; - } + if (m_ImageNode.IsNull()) + { + QMessageBox::information( nullptr, "Warning", "Please load and select a qball image before starting image processing."); + return; + } - if (dynamic_cast(m_ImageNode->GetData())) - m_QBallImage = dynamic_cast(m_ImageNode->GetData()); - else if (dynamic_cast(m_ImageNode->GetData())) - m_TensorImage = dynamic_cast(m_ImageNode->GetData()); + if (dynamic_cast(m_ImageNode->GetData())) + m_QBallImage = dynamic_cast(m_ImageNode->GetData()); + else if (dynamic_cast(m_ImageNode->GetData())) + m_TensorImage = dynamic_cast(m_ImageNode->GetData()); - if (m_QBallImage.IsNull() && m_TensorImage.IsNull()) - return; + if (m_QBallImage.IsNull() && m_TensorImage.IsNull()) + return; - // cast qbi to itk - m_TrackingNode = m_ImageNode; - m_ItkTensorImage = nullptr; - m_ItkQBallImage = nullptr; - m_MaskImage = nullptr; + // cast qbi to itk + m_TrackingNode = m_ImageNode; + m_ItkTensorImage = nullptr; + m_ItkQBallImage = nullptr; + m_MaskImage = nullptr; - if (m_QBallImage.IsNotNull()) - { - m_ItkQBallImage = ItkQBallImgType::New(); - mitk::CastToItkImage(m_QBallImage, m_ItkQBallImage); - } - else - { - m_ItkTensorImage = ItkTensorImage::New(); - mitk::CastToItkImage(m_TensorImage, m_ItkTensorImage); - } + if (m_QBallImage.IsNotNull()) + { + m_ItkQBallImage = ItkQBallImgType::New(); + mitk::CastToItkImage(m_QBallImage, m_ItkQBallImage); + } + else + { + m_ItkTensorImage = ItkTensorImage::New(); + mitk::CastToItkImage(m_TensorImage, m_ItkTensorImage); + } - // mask image found? - // catch exceptions thrown by the itkAccess macros - try{ - if(m_MaskImageNode.IsNotNull()) - { - if (dynamic_cast(m_MaskImageNode->GetData())) - mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); - } + // mask image found? + // catch exceptions thrown by the itkAccess macros + try{ + if(m_MaskImageNode.IsNotNull()) + { + if (dynamic_cast(m_MaskImageNode->GetData())) + mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); } - catch(...){}; + } + catch(...){}; - // start worker thread - m_TrackingThread.start(QThread::LowestPriority); + // start worker thread + m_TrackingThread.start(QThread::LowestPriority); } // generate mitkFiberBundle from tracking filter output void QmitkGibbsTrackingView::GenerateFiberBundle() { - if (m_GlobalTracker.IsNull() || (!(m_Controls->m_VisualizationCheckbox->isChecked() || m_Controls->m_VisualizeOnceButton->isChecked()) && m_ThreadIsRunning)) - return; + if (m_GlobalTracker.IsNull() || (!(m_Controls->m_VisualizationCheckbox->isChecked() || m_Controls->m_VisualizeOnceButton->isChecked()) && m_ThreadIsRunning)) + return; - if (m_Controls->m_VisualizeOnceButton->isChecked()) - m_Controls->m_VisualizeOnceButton->setChecked(false); + if (m_Controls->m_VisualizeOnceButton->isChecked()) + m_Controls->m_VisualizeOnceButton->setChecked(false); - vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); - if ( m_GlobalTracker->GetNumAcceptedFibers()==0 ) - return; - m_FiberBundle = mitk::FiberBundle::New(fiberBundle); - m_FiberBundle->SetReferenceGeometry(dynamic_cast(m_ImageNode->GetData())->GetGeometry()); + vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); + if ( m_GlobalTracker->GetNumAcceptedFibers()==0 ) + return; + m_FiberBundle = mitk::FiberBundle::New(fiberBundle); + m_FiberBundle->SetReferenceGeometry(dynamic_cast(m_ImageNode->GetData())->GetGeometry()); - if (m_FiberBundleNode.IsNotNull()){ - GetDataStorage()->Remove(m_FiberBundleNode); - m_FiberBundleNode = 0; - } - m_FiberBundleNode = mitk::DataNode::New(); - m_FiberBundleNode->SetData(m_FiberBundle); + if (m_FiberBundleNode.IsNotNull()){ + GetDataStorage()->Remove(m_FiberBundleNode); + m_FiberBundleNode = 0; + } + m_FiberBundleNode = mitk::DataNode::New(); + m_FiberBundleNode->SetData(m_FiberBundle); - QString name("FiberBundle_"); - name += m_ImageNode->GetName().c_str(); - name += "_Gibbs"; - m_FiberBundleNode->SetName(name.toStdString()); - m_FiberBundleNode->SetVisibility(true); + QString name("FiberBundle_"); + name += m_ImageNode->GetName().c_str(); + name += "_Gibbs"; + m_FiberBundleNode->SetName(name.toStdString()); + m_FiberBundleNode->SetVisibility(true); - if (!m_OutputFileName.isEmpty() && !m_ThreadIsRunning) + if (!m_OutputFileName.isEmpty() && !m_ThreadIsRunning) + { + try { - try - { - mitk::IOUtil::Save(m_FiberBundle.GetPointer(),m_OutputFileName.toStdString()); - QMessageBox::information(nullptr, "Fiber bundle saved to", m_OutputFileName); - } - catch (itk::ExceptionObject &ex) - { - QMessageBox::information(nullptr, "Fiber bundle could not be saved", QString("%1\n%2\n%3\n%4\n%5\n%6").arg(ex.GetNameOfClass()).arg(ex.GetFile()).arg(ex.GetLine()).arg(ex.GetLocation()).arg(ex.what()).arg(ex.GetDescription())); - } + mitk::IOUtil::Save(m_FiberBundle.GetPointer(),m_OutputFileName.toStdString()); + QMessageBox::information(nullptr, "Fiber bundle saved to", m_OutputFileName); } - if(m_ImageNode.IsNull()) - GetDataStorage()->Add(m_FiberBundleNode); - else - GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); + catch (itk::ExceptionObject &ex) + { + QMessageBox::information(nullptr, "Fiber bundle could not be saved", QString("%1\n%2\n%3\n%4\n%5\n%6").arg(ex.GetNameOfClass()).arg(ex.GetFile()).arg(ex.GetLine()).arg(ex.GetLocation()).arg(ex.what()).arg(ex.GetDescription())); + } + } + if(m_ImageNode.IsNull()) + GetDataStorage()->Add(m_FiberBundleNode); + else + GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); } void QmitkGibbsTrackingView::SetOutputFile() { - // SELECT FOLDER DIALOG - m_OutputFileName = QFileDialog::getSaveFileName(0, - tr("Set file name"), - QDir::currentPath()+"/FiberBundle.fib", - tr("Fiber Bundle (*.fib)") ); - if (m_OutputFileName.isEmpty()) - m_Controls->m_OutputFileLabel->setText("N/A"); - else - m_Controls->m_OutputFileLabel->setText(m_OutputFileName); + // SELECT FOLDER DIALOG + m_OutputFileName = QFileDialog::getSaveFileName(0, + tr("Set file name"), + QDir::currentPath()+"/FiberBundle.fib", + tr("Fiber Bundle (*.fib)") ); + if (m_OutputFileName.isEmpty()) + m_Controls->m_OutputFileLabel->setText("N/A"); + else + m_Controls->m_OutputFileLabel->setText(m_OutputFileName); } // save current tracking paramters as xml file (.gtp) void QmitkGibbsTrackingView::SaveTrackingParameters() { - TiXmlDocument documentXML; - TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); - documentXML.LinkEndChild( declXML ); - - TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); - mainXML->SetAttribute("file_version", "0.1"); - documentXML.LinkEndChild(mainXML); - - TiXmlElement* paramXML = new TiXmlElement("parameter_set"); - paramXML->SetAttribute("iterations", m_Controls->m_IterationsBox->text().toStdString()); - paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); - paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); - paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); - paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); - paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); - paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); - paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); - paramXML->SetAttribute("curvature_threshold", QString::number(m_Controls->m_CurvatureThresholdSlider->value()).toStdString()); - mainXML->LinkEndChild(paramXML); - QString filename = QFileDialog::getSaveFileName( - 0, - tr("Save Parameters"), - QDir::currentPath()+"/param.gtp", - tr("Global Tracking Parameters (*.gtp)") ); - - if(filename.isEmpty() || filename.isNull()) - return; - if(!filename.endsWith(".gtp")) - filename += ".gtp"; - documentXML.SaveFile( filename.toStdString() ); + TiXmlDocument documentXML; + TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); + documentXML.LinkEndChild( declXML ); + + TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); + mainXML->SetAttribute("file_version", "0.1"); + documentXML.LinkEndChild(mainXML); + + TiXmlElement* paramXML = new TiXmlElement("parameter_set"); + paramXML->SetAttribute("iterations", m_Controls->m_IterationsBox->text().toStdString()); + paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); + paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); + paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); + paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); + paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); + paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); + paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); + paramXML->SetAttribute("curvature_threshold", QString::number(m_Controls->m_CurvatureThresholdSlider->value()).toStdString()); + mainXML->LinkEndChild(paramXML); + QString filename = QFileDialog::getSaveFileName( + 0, + tr("Save Parameters"), + QDir::currentPath()+"/param.gtp", + tr("Global Tracking Parameters (*.gtp)") ); + + if(filename.isEmpty() || filename.isNull()) + return; + if(!filename.endsWith(".gtp")) + filename += ".gtp"; + documentXML.SaveFile( filename.toStdString() ); } // load current tracking paramters from xml file (.gtp) void QmitkGibbsTrackingView::LoadTrackingParameters() { - QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); - if(filename.isEmpty() || filename.isNull()) - return; - - TiXmlDocument doc( filename.toStdString() ); - doc.LoadFile(); - - TiXmlHandle hDoc(&doc); - TiXmlElement* pElem; - TiXmlHandle hRoot(0); - - pElem = hDoc.FirstChildElement().Element(); - hRoot = TiXmlHandle(pElem); - pElem = hRoot.FirstChildElement("parameter_set").Element(); - - QString iterations(pElem->Attribute("iterations")); - m_Controls->m_IterationsBox->setText(iterations); - - QString particleLength(pElem->Attribute("particle_length")); - float pLength = particleLength.toFloat(); - QString particleWidth(pElem->Attribute("particle_width")); - float pWidth = particleWidth.toFloat(); - - if (pLength==0) - m_Controls->m_ParticleLengthLabel->setText("auto"); - else - m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); - if (pWidth==0) - m_Controls->m_ParticleWidthLabel->setText("auto"); - else - m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); - - m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); - m_Controls->m_ParticleLengthSlider->setValue(pLength*10); - - QString partWeight(pElem->Attribute("particle_weight")); - m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); - m_Controls->m_ParticleWeightLabel->setText(partWeight); - - QString startTemp(pElem->Attribute("temp_start")); - m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); - m_Controls->m_StartTempLabel->setText(startTemp); - - QString endTemp(pElem->Attribute("temp_end")); - m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); - m_Controls->m_EndTempLabel->setText(endTemp); - - QString inExBalance(pElem->Attribute("inexbalance")); - m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); - m_Controls->m_InExBalanceLabel->setText(inExBalance); - - QString fiberLength(pElem->Attribute("fiber_length")); - m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); - m_Controls->m_FiberLengthLabel->setText(fiberLength+"mm"); - - QString curvThres(pElem->Attribute("curvature_threshold")); - m_Controls->m_CurvatureThresholdSlider->setValue(curvThres.toInt()); - m_Controls->m_CurvatureThresholdLabel->setText(curvThres+"°"); + QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); + if(filename.isEmpty() || filename.isNull()) + return; + + TiXmlDocument doc( filename.toStdString() ); + doc.LoadFile(); + + TiXmlHandle hDoc(&doc); + TiXmlElement* pElem; + TiXmlHandle hRoot(0); + + pElem = hDoc.FirstChildElement().Element(); + hRoot = TiXmlHandle(pElem); + pElem = hRoot.FirstChildElement("parameter_set").Element(); + + QString iterations(pElem->Attribute("iterations")); + m_Controls->m_IterationsBox->setText(iterations); + + QString particleLength(pElem->Attribute("particle_length")); + float pLength = particleLength.toFloat(); + QString particleWidth(pElem->Attribute("particle_width")); + float pWidth = particleWidth.toFloat(); + + if (pLength==0) + m_Controls->m_ParticleLengthLabel->setText("auto"); + else + m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); + if (pWidth==0) + m_Controls->m_ParticleWidthLabel->setText("auto"); + else + m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); + + m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); + m_Controls->m_ParticleLengthSlider->setValue(pLength*10); + + QString partWeight(pElem->Attribute("particle_weight")); + m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); + m_Controls->m_ParticleWeightLabel->setText(partWeight); + + QString startTemp(pElem->Attribute("temp_start")); + m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); + m_Controls->m_StartTempLabel->setText(startTemp); + + QString endTemp(pElem->Attribute("temp_end")); + m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); + m_Controls->m_EndTempLabel->setText(endTemp); + + QString inExBalance(pElem->Attribute("inexbalance")); + m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); + m_Controls->m_InExBalanceLabel->setText(inExBalance); + + QString fiberLength(pElem->Attribute("fiber_length")); + m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); + m_Controls->m_FiberLengthLabel->setText(fiberLength+"mm"); + + QString curvThres(pElem->Attribute("curvature_threshold")); + m_Controls->m_CurvatureThresholdSlider->setValue(curvThres.toInt()); + m_Controls->m_CurvatureThresholdLabel->setText(curvThres+"°"); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingViewControls.ui index d43fb6ad0f..ce2353faad 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkGibbsTrackingViewControls.ui @@ -1,1103 +1,1121 @@ QmitkGibbsTrackingViewControls 0 0 463 1011 0 0 0 0 QmitkTemplate QFormLayout::AllNonFixedFieldsGrow - - - - Please Select Input Data - - - - - - Q-Ball/Tensor Image: - - - - - - - Mandatory input - - - <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> - - - true - - - - - - - Mask Image: - - - - - - - <html><head/><body><p>White matter probability mask image.</p></body></html> - - - <html><head/><body><p><span style=" color:#969696;">optional</span></p></body></html> - - - true - - - - - - QFrame::NoFrame QFrame::Plain 0 0 0 0 0 0 false No Q-Ball image selected. Qt::LeftToRight Start Tractography :/QtWidgetsExt/play.xpm:/QtWidgetsExt/play.xpm false Qt::LeftToRight Stop Tractography :/QtWidgetsExt/stop.xpm:/QtWidgetsExt/stop.xpm Parameters 0 Output File: Advanced Settings true Activate continuous visualization of intermediate results. Visualize Tractography true QFrame::NoFrame QFrame::Plain 0 0 0 0 0 0 Select output file name and folder. ... N/A true Visualize intermediate result. - + :/QmitkDiffusionImaging/Refresh_48.png:/QmitkDiffusionImaging/Refresh_48.png true Iterations: 1e9 true QFrame::StyledPanel QFrame::Raised 9 0 9 0 4 Particle Width: 0 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 0.1 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter Particle Weight: 1 99 1 10 Qt::Horizontal QSlider::NoTicks Start Temperature: automatic estimation from gfa map and q-ball data. 0 1000 1 0 Qt::Horizontal true QSlider::NoTicks IE Bias < 0 < EE Bias -50 50 1 Qt::Horizontal QSlider::NoTicks 0.001 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter Curvature Threshold: Balance In/Ex Energy: 45° Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto = 0.5 * min. spacing; sigma 100 1 Qt::Horizontal QSlider::NoTicks Particle Length: auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter Min. Fiber Length: Only fibers longer than specified are accepted. 500 1 20 Qt::Horizontal QSlider::NoTicks Allow only fiber curvature values smaller than the selected threshold. 180 1 45 Qt::Horizontal QSlider::NoTicks End Temperature: 20mm Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 1 100 1 10 Qt::Horizontal false false QSlider::NoTicks auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto = 1.5 * min. spacing; l 100 1 Qt::Horizontal QSlider::NoTicks Random Seed auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto = 1.5 * min. spacing; l -1 100 1 -1 Qt::Horizontal QSlider::NoTicks QFrame::NoFrame QFrame::Plain 0 0 0 0 0 0 true Save current parameters as xml (.gtp) Qt::LeftToRight Save Parameters :/QtWidgetsExt/btnMoveDown.png:/QtWidgetsExt/btnMoveDown.png true Load parameters from xml file (.gtp) Qt::LeftToRight Load Parameters :/QtWidgetsExt/btnMoveUp.png:/QtWidgetsExt/btnMoveUp.png - + + + + Qt::Vertical + + + QSizePolicy::Expanding + + + + 0 + 0 + + + + + + + + + 0 + 0 + + + + Please Select Input Data + + + + + + + 0 + 0 + + + + Q-Ball/Tensor Image: + + + + + + + Mandatory input + + + <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> + + + true + + + + + + + Mask Image: + + + + + + + <html><head/><body><p>White matter probability mask image.</p></body></html> + + + <html><head/><body><p><span style=" color:#969696;">optional</span></p></body></html> + + + true + + + + + + + + + + 0 + 0 + + Monitor Progress: - Will only be updated if tracking is visualized Will only be updated if tracking is visualized Accepted Fibers: Connections: Particles: Proposal Acceptance Rate: Tracking Time: Will only be updated if tracking is visualized - - - - - - - - - Qt::Vertical - - - QSizePolicy::Expanding - - - - 0 - 0 - - - - - + 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 d83356dee7..503be0b386 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,747 +1,819 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include #include #include #include #include #include #include #include #include #include #include // VTK #include #include #include #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; +QmitkStreamlineTrackingWorker::QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view) + : m_View(view) +{ +} + +void QmitkStreamlineTrackingWorker::run() +{ + m_View->m_Tracker->Update(); + m_View->m_TrackingThread.quit(); +} + QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() - : m_Controls(nullptr) + : m_TrackingWorker(this) + , m_Controls(nullptr) , m_TrackingHandler(nullptr) + , m_ThreadIsRunning(false) { + m_TrackingWorker.moveToThread(&m_TrackingThread); + connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); + connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); + connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); + m_TrackingTimer = new QTimer(this); } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { + if (m_Tracker.IsNull()) + return; + m_Tracker->SetStopTracking(true); + m_TrackingThread.wait(); } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TissueImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ForestBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestBox->SetPredicate(isTractographyForest); m_Controls->m_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageBox->SetZeroEntryText("--"); m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_SeedImageBox->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_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); + connect( m_Controls->commandLinkButton_2, SIGNAL(clicked()), this, SLOT(StopTractography()) ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_TissueImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(valueChanged(double)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(valueChanged(int)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedGmBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); + m_FirstTensorProbRun = true; + StartStopTrackingGui(false); + } + UpdateGui(); +} - m_FirstTensorProbRun = true; +void QmitkStreamlineTrackingView::StopTractography() +{ + if (m_Tracker.IsNull()) + return; + + m_Tracker->SetStopTracking(true); +} + +void QmitkStreamlineTrackingView::TimerUpdate() +{ + if (m_Tracker.IsNull()) + return; + + QString status_text(m_Tracker->GetStatusText().c_str()); + m_Controls->m_StatusTextBox->setText(status_text); +} + +void QmitkStreamlineTrackingView::BeforeThread() +{ + m_TrackingTimer->start(1000); +} + +void QmitkStreamlineTrackingView::AfterThread() +{ + m_TrackingTimer->stop(); + if (!m_Tracker->GetUseOutputProbabilityMap()) + { + vtkSmartPointer fiberBundle = m_Tracker->GetFiberPolyData(); + + if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) + { + QMessageBox warnBox; + warnBox.setWindowTitle("Warning"); + warnBox.setText("No fiberbundle was generated!"); + warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); + warnBox.setIcon(QMessageBox::Warning); + warnBox.exec(); + 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(); + + if (m_Controls->m_InteractiveBox->isChecked()) + { + if (m_InteractiveNode.IsNull()) + { + m_InteractiveNode = mitk::DataNode::New(); + QString name("Interactive"); + m_InteractiveNode->SetName(name.toStdString()); + GetDataStorage()->Add(m_InteractiveNode); + } + m_InteractiveNode->SetData(fib); + + if (auto renderWindowPart = this->GetRenderWindowPart()) + renderWindowPart->RequestUpdate(); + } + else + { + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(fib); + QString name("FiberBundle_"); + name += m_InputImageNodes.at(0)->GetName().c_str(); + name += "_Streamline"; + node->SetName(name.toStdString()); + GetDataStorage()->Add(node, m_InputImageNodes.at(0)); + } } + else + { + TrackerType::ItkDoubleImgType::Pointer outImg = m_Tracker->GetOutputProbabilityMap(); + mitk::Image::Pointer img = mitk::Image::New(); + img->InitializeByItk(outImg.GetPointer()); + img->SetVolume(outImg->GetBufferPointer()); + + if (m_Controls->m_InteractiveBox->isChecked()) + { + if (m_InteractiveNode.IsNull()) + { + m_InteractiveNode = mitk::DataNode::New(); + QString name("Interactive"); + m_InteractiveNode->SetName(name.toStdString()); + GetDataStorage()->Add(m_InteractiveNode); + } + m_InteractiveNode->SetData(img); + + mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); + lut->SetType(mitk::LookupTable::JET_TRANSPARENT); + mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); + lut_prop->SetLookupTable(lut); + m_InteractiveNode->SetProperty("LookupTable", lut_prop); + m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); + + if (auto renderWindowPart = this->GetRenderWindowPart()) + renderWindowPart->RequestUpdate(); + } + else + { + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(img); + QString name("ProbabilityMap_"); + name += m_InputImageNodes.at(0)->GetName().c_str(); + node->SetName(name.toStdString()); + + mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); + lut->SetType(mitk::LookupTable::JET_TRANSPARENT); + mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); + lut_prop->SetLookupTable(lut); + node->SetProperty("LookupTable", lut_prop); + node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); + GetDataStorage()->Add(node, m_InputImageNodes.at(0)); + } + } + StartStopTrackingGui(false); UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked())) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } DoFiberTracking(); } void QmitkStreamlineTrackingView::OnParameterChanged() { if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); mitk::IRenderWindowPart* renderWindow = this->GetRenderWindowPart(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedGmBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_TissueImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_10->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); if (renderWindow) { { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("axial"))->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); m_SliceObserverTag1 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); } { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("sagittal"))->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); m_SliceObserverTag2 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); } { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("coronal"))->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); m_SliceObserverTag3 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); } } } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; if (renderWindow) { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("axial"))->GetSliceNavigationController(); slicer->RemoveObserver(m_SliceObserverTag1); slicer = renderWindow->GetQmitkRenderWindow(QString("sagittal"))->GetSliceNavigationController(); slicer->RemoveObserver(m_SliceObserverTag2); slicer = renderWindow->GetQmitkRenderWindow(QString("coronal"))->GetSliceNavigationController(); slicer->RemoveObserver(m_SliceObserverTag3); } } } void QmitkStreamlineTrackingView::OnSliceChanged(const itk::EventObject& /*e*/) { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (m_TrackingHandler != nullptr) { delete m_TrackingHandler; m_TrackingHandler = nullptr; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { m_InputImageNodes.clear(); m_InputImages.clear(); m_AdditionalInputImages.clear(); DeleteTrackingHandler(); 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 if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( 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())); } else if (dim==3) { m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } } } UpdateGui(); OnParameterChanged(); } void QmitkStreamlineTrackingView::UpdateGui() { m_Controls->m_TensorImageLabel->setText("mandatory"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageBox->setEnabled(false); m_Controls->mFaImageLabel->setEnabled(false); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestBox->setEnabled(false); m_Controls->m_ForestLabel->setEnabled(false); + m_Controls->commandLinkButton->setEnabled(false); if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_SeedGmBox->setEnabled(true); else m_Controls->m_SeedGmBox->setEnabled(false); if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText(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(!m_Controls->m_InteractiveBox->isChecked()); + m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked() && !m_ThreadIsRunning); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestBox->setEnabled(true); m_Controls->m_ForestLabel->setEnabled(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } else - { m_Controls->m_InputData->setTitle("Please Select Input Data"); - m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); - } } +void QmitkStreamlineTrackingView::StartStopTrackingGui(bool start) +{ + m_ThreadIsRunning = start; + + if (!m_Controls->m_InteractiveBox->isChecked()) + { + m_Controls->commandLinkButton_2->setVisible(start); + m_Controls->commandLinkButton->setVisible(!start); + m_Controls->m_InteractiveBox->setEnabled(!start); + m_Controls->m_StatusTextBox->setVisible(start); + } +} + void QmitkStreamlineTrackingView::DoFiberTracking() { + if (m_ThreadIsRunning) + return; if (m_InputImages.empty()) return; + if (m_Controls->m_InteractiveBox->isChecked() && m_SeedPoints.empty()) + return; + StartStopTrackingGui(true); - typedef itk::StreamlineTrackingFilter TrackerType; - TrackerType::Pointer tracker = TrackerType::New(); + m_Tracker = TrackerType::New(); 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!"); + StartStopTrackingGui(false); return; } if (m_FirstTensorProbRun) { QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. ODFs are sharpened (min-max normalized and raised to the power of 4). TEND parameters are ignored."); m_FirstTensorProbRun = false; } if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); 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(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); } else { if (m_TrackingHandler==nullptr) { m_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(m_TrackingHandler)->AddTensorImage(itkImg); } if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(m_TrackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); dynamic_cast(m_TrackingHandler)->SetOdfImage(itkImg); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestBox->GetSelectedNode().IsNull() ) { QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography selected."); + StartStopTrackingGui(false); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestBox->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { int num_previous_directions = (forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { int num_previous_directions = (forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); + StartStopTrackingGui(false); return; } } } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); + StartStopTrackingGui(false); return; } try { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(m_InputImages.at(0)); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); m_TrackingHandler = new mitk::TrackingHandlerPeaks(); dynamic_cast(m_TrackingHandler)->SetPeakImage(itkImg); } dynamic_cast(m_TrackingHandler)->SetPeakThreshold(m_Controls->m_ScalarThresholdBox->value()); } catch(...) { + StartStopTrackingGui(false); return; } } m_TrackingHandler->SetFlipX(m_Controls->m_FlipXBox->isChecked()); m_TrackingHandler->SetFlipY(m_Controls->m_FlipYBox->isChecked()); m_TrackingHandler->SetFlipZ(m_Controls->m_FlipZBox->isChecked()); m_TrackingHandler->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); switch (m_Controls->m_ModeBox->currentIndex()) { case 0: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); break; case 1: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); break; default: m_TrackingHandler->SetMode(mitk::TrackingDataHandler::MODE::DETERMINISTIC); } if (m_Controls->m_InteractiveBox->isChecked()) { - tracker->SetSeedPoints(m_SeedPoints); + m_Tracker->SetSeedPoints(m_SeedPoints); } else 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); + m_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); + m_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); + m_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(m_Controls->m_SamplingDistanceBox->value()); - tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); - tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); - tracker->SetAposterioriCurvCheck(false); - tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); - tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); - tracker->SetTrackingHandler(m_TrackingHandler); - tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); - tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); - tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); - tracker->Update(); - - if (!tracker->GetUseOutputProbabilityMap()) - { - vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); - - if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) - { - QMessageBox warnBox; - warnBox.setWindowTitle("Warning"); - warnBox.setText("No fiberbundle was generated!"); - warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); - warnBox.setIcon(QMessageBox::Warning); - warnBox.exec(); - 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(); - - if (m_Controls->m_InteractiveBox->isChecked()) - { - if (m_InteractiveNode.IsNull()) - { - m_InteractiveNode = mitk::DataNode::New(); - QString name("Interactive"); - m_InteractiveNode->SetName(name.toStdString()); - GetDataStorage()->Add(m_InteractiveNode); - } - m_InteractiveNode->SetData(fib); - - if (auto renderWindowPart = this->GetRenderWindowPart()) - renderWindowPart->RequestUpdate(); - } - else - { - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(fib); - QString name("FiberBundle_"); - name += m_InputImageNodes.at(0)->GetName().c_str(); - name += "_Streamline"; - node->SetName(name.toStdString()); - GetDataStorage()->Add(node, m_InputImageNodes.at(0)); - } + m_Tracker->SetTissueImage(mask); + m_Tracker->SetSeedOnlyGm(m_Controls->m_SeedGmBox->isChecked()); } - else - { - TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); - mitk::Image::Pointer img = mitk::Image::New(); - img->InitializeByItk(outImg.GetPointer()); - img->SetVolume(outImg->GetBufferPointer()); - if (m_Controls->m_InteractiveBox->isChecked()) - { - if (m_InteractiveNode.IsNull()) - { - m_InteractiveNode = mitk::DataNode::New(); - QString name("Interactive"); - m_InteractiveNode->SetName(name.toStdString()); - GetDataStorage()->Add(m_InteractiveNode); - } - m_InteractiveNode->SetData(img); - - mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); - lut->SetType(mitk::LookupTable::JET_TRANSPARENT); - mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); - lut_prop->SetLookupTable(lut); - m_InteractiveNode->SetProperty("LookupTable", lut_prop); - m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); - - if (auto renderWindowPart = this->GetRenderWindowPart()) - renderWindowPart->RequestUpdate(); - } - else - { - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(img); - QString name("ProbabilityMap_"); - name += m_InputImageNodes.at(0)->GetName().c_str(); - node->SetName(name.toStdString()); - - mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); - lut->SetType(mitk::LookupTable::JET_TRANSPARENT); - mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); - lut_prop->SetLookupTable(lut); - node->SetProperty("LookupTable", lut_prop); - node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); - - GetDataStorage()->Add(node, m_InputImageNodes.at(0)); - } - } + m_Tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); + m_Tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); + m_Tracker->SetSamplingDistance(m_Controls->m_SamplingDistanceBox->value()); + m_Tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); + m_Tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); + m_Tracker->SetAposterioriCurvCheck(false); + m_Tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); + m_Tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); + m_Tracker->SetTrackingHandler(m_TrackingHandler); + m_Tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); + m_Tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); + m_Tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); + + m_TrackingThread.start(QThread::LowestPriority); } 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 a2a2d98fc8..4c377537d2 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,107 +1,140 @@ /*=================================================================== 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 #include #include #include #include #include #include +#include +#include +class QmitkStreamlineTrackingView; + +class QmitkStreamlineTrackingWorker : public QObject +{ + Q_OBJECT + +public: + + QmitkStreamlineTrackingWorker(QmitkStreamlineTrackingView* view); + +public slots: + + void run(); + +private: + + QmitkStreamlineTrackingView* m_View; +}; /*! \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; + typedef itk::Image< unsigned char, 3 > ItkUCharImageType; + typedef itk::Image< float, 3 > ItkFloatImageType; + typedef itk::StreamlineTrackingFilter TrackerType; QmitkStreamlineTrackingView(); virtual ~QmitkStreamlineTrackingView(); virtual void CreateQtPartControl(QWidget *parent) override; /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; + TrackerType::Pointer m_Tracker; + QmitkStreamlineTrackingWorker m_TrackingWorker; + QThread m_TrackingThread; + protected slots: void DoFiberTracking(); ///< start fiber tracking void UpdateGui(); void ToggleInteractive(); void DeleteTrackingHandler(); void OnParameterChanged(); void InteractiveSeedChanged(bool posChanged=false); void ForestSwitched(); void OutputStyleSwitched(); + void AfterThread(); ///< update gui etc. after tracking has finished + void BeforeThread(); ///< start timer etc. + void TimerUpdate(); + void StopTractography(); 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: + void StartStopTrackingGui(bool start); + void OnSliceChanged(const itk::EventObject& e); + int m_SliceObserverTag1; int m_SliceObserverTag2; int m_SliceObserverTag3; std::vector< itk::Point > m_SeedPoints; mitk::DataNode::Pointer m_InteractiveNode; mitk::DataNode::Pointer m_InteractivePointSetNode; std::vector< mitk::DataNode::Pointer > m_InputImageNodes; ///< input image nodes std::vector< mitk::Image::Pointer > m_InputImages; ///< input images std::vector< mitk::Image::Pointer > m_AdditionalInputImages; bool m_FirstTensorProbRun; - void OnSliceChanged(const itk::EventObject& e); - mitk::TrackingDataHandler* m_TrackingHandler; + mitk::TrackingDataHandler* m_TrackingHandler; + bool m_ThreadIsRunning; + QTimer* m_TrackingTimer; }; #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 4304dcc8ac..17129d7efd 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,1033 +1,1059 @@ QmitkStreamlineTrackingViewControls 0 0 413 1222 0 0 QmitkTemplate 0 0 3 3 Please Select Input Data Mask Image: Fibers that enter a region defined in this image will stop immediately. QComboBox::AdjustToMinimumContentsLength - Seed points are only placed inside the mask image. If no seed mask is selected, the whole image is seeded. QComboBox::AdjustToMinimumContentsLength - Stop 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 - Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. QComboBox::AdjustToMinimumContentsLength - FA/GFA Image: If an image is selected, the stopping criterion is not calculated from the input image but instead the selected image is used. QComboBox::AdjustToMinimumContentsLength - Seed Image: Tissue Image: Input Image. ODF, tensor, peak, and, in case of ML tractography, raw diffusion-weighted images are currently supported. <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Input Image. ODF, tensor and peak images are currently supported. Input Image: Tractography Forest: Random forest for machine learning based tractography. QComboBox::AdjustToMinimumContentsLength - - + false Start Tractography - + 0 0 Interactive Tractography Number of seed points normally distributed around selected position. 1 9999999 50 Num. Seeds: true Dynamically pick seed location by click into image. Enable Interactive Tractography Seedpoints are normally distributed within a sphere centered at the selected position with the specified radius (in mm). 2 50.000000000000000 0.100000000000000 2.500000000000000 Radius: true When checked, parameter changes cause instant retracking while in interactive mode. Update on Parameter Change true - + 0 0 Parameters Mode: Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. Deterministic Probabilistic Seeds per Voxel: Number of seed points placed in each voxel. 1 9999999 Max. num. fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. -1 999999999 -1 Cutoff: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 ODF Cutoff: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. 5 1.000000000000000 0.100000000000000 0.100000000000000 If you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). This is not necessary for CSD fODFs, since they are naturally much sharper. Sharpen ODFs Qt::Horizontal QSizePolicy::Fixed 200 0 - + Advanced Parameters QFrame::NoFrame QFrame::Raised 0 0 0 0 Min. Tract Length: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Angular Threshold: Flip directions: f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 1.000000000000000 f parameter of tensor tractography. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). f: g: Shorter fibers are discarded. Minimum fiber length (in mm) 1 999.000000000000000 1.000000000000000 20.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 Internally flips progression directions. This might be necessary depending on the input data. x Internally flips progression directions. This might be necessary depending on the input data. y Internally flips progression directions. This might be necessary depending on the input data. z Step Size: Default: 90° * step_size -1 90 1 -1 f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 0.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 Requires tissue image. Enable Gray Matter Seeding false If false, nearest neighbor interpolation is used. Enable Trilinear Interpolation true - + 0 0 Neighborhood Sampling QFrame::NoFrame QFrame::Raised 0 0 0 0 Num. Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 Sampling Distance: Sampling distance (in voxels) 2 10.000000000000000 0.100000000000000 0.250000000000000 Only neighborhood samples in front of the current streamline position are considered. Use Only Frontal Samples false If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. Use Stop-Votes false - + Qt::Vertical QSizePolicy::Expanding 20 220 - + 0 0 Output and Postprocessing QFrame::NoFrame QFrame::Raised 0 0 0 0 Resample fibers using the specified error constraint. Compress Fibers true Qt::StrongFocus Lossy fiber compression. Recommended for large tractograms. Maximum error in mm. 3 10.000000000000000 0.010000000000000 0.100000000000000 Output probability map instead of tractogram. Output Probability Map false + + + + true + + + Stop Tractography + + + + + + + true + + + + 0 + 0 + + + + true + + + QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h