diff --git a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkNrrdTensorImageReader.cpp b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkNrrdTensorImageReader.cpp index 471d996649..52187ec588 100644 --- a/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkNrrdTensorImageReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/autoload/IO/mitkNrrdTensorImageReader.cpp @@ -1,420 +1,421 @@ /*=================================================================== 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 "mitkNrrdTensorImageReader.h" #include #include "mitkDiffusionCoreIOMimeTypes.h" #include "itkImageFileReader.h" #include "itkImageRegionIterator.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" #include #include "mitkITKImageImport.h" #include "mitkImageDataItem.h" #include +#include namespace mitk { NrrdTensorImageReader::NrrdTensorImageReader(const NrrdTensorImageReader& other) : mitk::AbstractFileReader(other) { } NrrdTensorImageReader::NrrdTensorImageReader() : mitk::AbstractFileReader( CustomMimeType( mitk::DiffusionCoreIOMimeTypes::DTI_MIMETYPE() ), mitk::DiffusionCoreIOMimeTypes::DTI_MIMETYPE_DESCRIPTION() ) { m_ServiceReg = this->RegisterService(); } NrrdTensorImageReader::~NrrdTensorImageReader() { } std::vector > NrrdTensorImageReader::Read() { std::vector > result; std::string location = GetInputLocation(); if ( location == "") { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, the filename is empty!"); } else { try { mitk::LocaleSwitch localeSwitch("C"); try { - std::string fname3 = "temp_dti.nii"; + std::string fname3 = mitk::IOUtil::GetTempPath()+"/temp_dti.nii"; itksys::SystemTools::CopyAFile(location.c_str(), fname3.c_str()); typedef itk::VectorImage ImageType; itk::NiftiImageIO::Pointer io = itk::NiftiImageIO::New(); typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetImageIO(io); reader->SetFileName(fname3); reader->Update(); ImageType::Pointer img = reader->GetOutput(); typedef itk::Image,3> VecImgType; VecImgType::Pointer vecImg = VecImgType::New(); vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin vecImg->SetDirection( img->GetDirection() ); // Set the image direction vecImg->SetRegions( img->GetLargestPossibleRegion()); vecImg->Allocate(); itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); ot.GoToBegin(); itk::ImageRegionIterator it (img, img->GetLargestPossibleRegion() ); it.GoToBegin(); typedef ImageType::PixelType VarPixType; typedef VecImgType::PixelType FixPixType; int numComponents = img->GetNumberOfComponentsPerPixel(); if (numComponents==6) { MITK_INFO << "Trying to load dti as 6-comp nifti ..."; while (!it.IsAtEnd()) { VarPixType vec = it.Get(); FixPixType fixVec(vec.GetDataPointer()); itk::DiffusionTensor3D tensor; tensor.SetElement(0, vec.GetElement(0)); tensor.SetElement(1, vec.GetElement(1)); tensor.SetElement(2, vec.GetElement(2)); tensor.SetElement(3, vec.GetElement(3)); tensor.SetElement(4, vec.GetElement(4)); tensor.SetElement(5, vec.GetElement(5)); fixVec = tensor; ot.Set(fixVec); ++ot; ++it; } } else if(numComponents==9) { MITK_INFO << "Trying to load dti as 9-comp nifti ..."; while (!it.IsAtEnd()) { VarPixType vec = it.Get(); itk::DiffusionTensor3D tensor; tensor.SetElement(0, vec.GetElement(0)); tensor.SetElement(1, vec.GetElement(1)); tensor.SetElement(2, vec.GetElement(2)); tensor.SetElement(3, vec.GetElement(4)); tensor.SetElement(4, vec.GetElement(5)); tensor.SetElement(5, vec.GetElement(8)); FixPixType fixVec(tensor); ot.Set(fixVec); ++ot; ++it; } } else if (numComponents==1) { MITK_INFO << "Trying to load dti as 4D nifti ..."; typedef itk::Image ImageType; typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetImageIO(io); reader->SetFileName(fname3); reader->Update(); ImageType::Pointer img = reader->GetOutput(); itk::Size<4> size = img->GetLargestPossibleRegion().GetSize(); while (!ot.IsAtEnd()) { itk::DiffusionTensor3D tensor; ImageType::IndexType idx; idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2]; if (size[3]==6) { for (unsigned int te=0; teGetPixel(idx)); } } else if (size[3]==9) { idx[3] = 0; tensor.SetElement(0, img->GetPixel(idx)); idx[3] = 1; tensor.SetElement(1, img->GetPixel(idx)); idx[3] = 2; tensor.SetElement(2, img->GetPixel(idx)); idx[3] = 4; tensor.SetElement(3, img->GetPixel(idx)); idx[3] = 5; tensor.SetElement(4, img->GetPixel(idx)); idx[3] = 8; tensor.SetElement(5, img->GetPixel(idx)); } else throw itk::ImageFileReaderException(__FILE__, __LINE__, "Unknown number of components for DTI file. Should be 6 or 9!"); FixPixType fixVec(tensor); ot.Set(fixVec); ++ot; } } OutputType::Pointer resultImage = OutputType::New(); resultImage->InitializeByItk( vecImg.GetPointer() ); resultImage->SetVolume( vecImg->GetBufferPointer() ); result.push_back( resultImage.GetPointer() ); } catch(...) { MITK_INFO << "Trying to load dti as nrrd ..."; typedef itk::VectorImage ImageType; itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetImageIO(io); reader->SetFileName(location); reader->Update(); ImageType::Pointer img = reader->GetOutput(); typedef itk::Image,3> VecImgType; VecImgType::Pointer vecImg = VecImgType::New(); vecImg->SetSpacing( img->GetSpacing() ); // Set the image spacing vecImg->SetOrigin( img->GetOrigin() ); // Set the image origin vecImg->SetDirection( img->GetDirection() ); // Set the image direction vecImg->SetRegions( img->GetLargestPossibleRegion()); vecImg->Allocate(); itk::ImageRegionIterator ot (vecImg, vecImg->GetLargestPossibleRegion() ); ot.GoToBegin(); itk::ImageRegionIterator it (img, img->GetLargestPossibleRegion() ); it.GoToBegin(); typedef ImageType::PixelType VarPixType; typedef VecImgType::PixelType FixPixType; int numComponents = img->GetNumberOfComponentsPerPixel(); itk::MetaDataDictionary imgMetaDictionary = img->GetMetaDataDictionary(); std::vector imgMetaKeys = imgMetaDictionary.GetKeys(); std::vector::const_iterator itKey = imgMetaKeys.begin(); std::string metaString; bool readFrame = false; double xx, xy, xz, yx, yy, yz, zx, zy, zz; MeasurementFrameType measFrame; measFrame.SetIdentity(); MeasurementFrameType measFrameTransp; measFrameTransp.SetIdentity(); for (; itKey != imgMetaKeys.end(); itKey ++) { itk::ExposeMetaData (imgMetaDictionary, *itKey, metaString); if (itKey->find("measurement frame") != std::string::npos) { sscanf(metaString.c_str(), " ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) ( %lf , %lf , %lf ) \n", &xx, &xy, &xz, &yx, &yy, &yz, &zx, &zy, &zz); if (xx>10e-10 || xy>10e-10 || xz>10e-10 || yx>10e-10 || yy>10e-10 || yz>10e-10 || zx>10e-10 || zy>10e-10 || zz>10e-10 ) { readFrame = true; measFrame(0,0) = xx; measFrame(0,1) = xy; measFrame(0,2) = xz; measFrame(1,0) = yx; measFrame(1,1) = yy; measFrame(1,2) = yz; measFrame(2,0) = zx; measFrame(2,1) = zy; measFrame(2,2) = zz; measFrameTransp = measFrame.GetTranspose(); } } } if (numComponents==6) { while (!it.IsAtEnd()) { // T'=RTR' VarPixType vec = it.Get(); FixPixType fixVec(vec.GetDataPointer()); if(readFrame) { itk::DiffusionTensor3D tensor; tensor.SetElement(0, vec.GetElement(0)); tensor.SetElement(1, vec.GetElement(1)); tensor.SetElement(2, vec.GetElement(2)); tensor.SetElement(3, vec.GetElement(3)); tensor.SetElement(4, vec.GetElement(4)); tensor.SetElement(5, vec.GetElement(5)); tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); fixVec = tensor; } ot.Set(fixVec); ++ot; ++it; } } else if(numComponents==9) { while (!it.IsAtEnd()) { VarPixType vec = it.Get(); itk::DiffusionTensor3D tensor; tensor.SetElement(0, vec.GetElement(0)); tensor.SetElement(1, vec.GetElement(1)); tensor.SetElement(2, vec.GetElement(2)); tensor.SetElement(3, vec.GetElement(4)); tensor.SetElement(4, vec.GetElement(5)); tensor.SetElement(5, vec.GetElement(8)); if(readFrame) { tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); } FixPixType fixVec(tensor); ot.Set(fixVec); ++ot; ++it; } } else if (numComponents==1) { typedef itk::Image ImageType; itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); typedef itk::ImageFileReader FileReaderType; FileReaderType::Pointer reader = FileReaderType::New(); reader->SetImageIO(io); reader->SetFileName(location); reader->Update(); ImageType::Pointer img = reader->GetOutput(); itk::Size<4> size = img->GetLargestPossibleRegion().GetSize(); while (!ot.IsAtEnd()) { itk::DiffusionTensor3D tensor; ImageType::IndexType idx; idx[0] = ot.GetIndex()[0]; idx[1] = ot.GetIndex()[1]; idx[2] = ot.GetIndex()[2]; if (size[3]==6) { for (unsigned int te=0; teGetPixel(idx)); } } else if (size[3]==9) { idx[3] = 0; tensor.SetElement(0, img->GetPixel(idx)); idx[3] = 1; tensor.SetElement(1, img->GetPixel(idx)); idx[3] = 2; tensor.SetElement(2, img->GetPixel(idx)); idx[3] = 4; tensor.SetElement(3, img->GetPixel(idx)); idx[3] = 5; tensor.SetElement(4, img->GetPixel(idx)); idx[3] = 8; tensor.SetElement(5, img->GetPixel(idx)); } else throw itk::ImageFileReaderException(__FILE__, __LINE__, "Unknown number of komponents for DTI file. Should be 6 or 9!"); if(readFrame) { tensor = ConvertMatrixTypeToFixedArrayType(tensor.PreMultiply(measFrame)); tensor = ConvertMatrixTypeToFixedArrayType(tensor.PostMultiply(measFrameTransp)); } FixPixType fixVec(tensor); ot.Set(fixVec); ++ot; } } else { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Image has wrong number of pixel components!"); } OutputType::Pointer resultImage = OutputType::New(); resultImage->InitializeByItk( vecImg.GetPointer() ); resultImage->SetVolume( vecImg->GetBufferPointer() ); result.push_back( resultImage.GetPointer() ); } } catch(std::exception& e) { throw itk::ImageFileReaderException(__FILE__, __LINE__, e.what()); } catch(...) { throw itk::ImageFileReaderException(__FILE__, __LINE__, "Sorry, an error occurred while reading the requested DTI file!"); } } return result; } itk::DiffusionTensor3D NrrdTensorImageReader ::ConvertMatrixTypeToFixedArrayType(const itk::DiffusionTensor3D::Superclass::MatrixType & matrix) { /* | 0 1 2 | * | X 3 4 | * | X X 5 | */ itk::DiffusionTensor3D arr; arr.SetElement(0,matrix(0,0)); arr.SetElement(1,matrix(0,1)); arr.SetElement(2,matrix(0,2)); arr.SetElement(3,matrix(1,3)); arr.SetElement(4,matrix(1,4)); arr.SetElement(5,matrix(2,5)); return arr; } } //namespace MITK mitk::NrrdTensorImageReader* mitk::NrrdTensorImageReader::Clone() const { return new NrrdTensorImageReader(*this); } diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkShCoefficientImageImporter.cpp b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkShCoefficientImageImporter.cpp index 6d2f924976..186146deb3 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkShCoefficientImageImporter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkShCoefficientImageImporter.cpp @@ -1,192 +1,192 @@ /*=================================================================== 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 __itkShCoefficientImageImporter_cpp #define __itkShCoefficientImageImporter_cpp #include #include #include #include "itkShCoefficientImageImporter.h" #include #include using namespace boost::math; namespace itk { template< class PixelType, int ShOrder > ShCoefficientImageImporter< PixelType, ShOrder >::ShCoefficientImageImporter() : m_Toolkit(FSL) { m_ShBasis.set_size(QBALL_ODFSIZE, (ShOrder+1)*(ShOrder+2)/2); } template< class PixelType, int ShOrder > void ShCoefficientImageImporter< PixelType, ShOrder > ::GenerateData() { CalcShBasis(); if (m_InputImage.IsNull()) return; Vector spacing4 = m_InputImage->GetSpacing(); Point origin4 = m_InputImage->GetOrigin(); Matrix direction4 = m_InputImage->GetDirection(); ImageRegion<4> imageRegion4 = m_InputImage->GetLargestPossibleRegion(); Vector spacing3; Point origin3; Matrix direction3; ImageRegion<3> imageRegion3; 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]; imageRegion3.SetSize(0, imageRegion4.GetSize()[0]); imageRegion3.SetSize(1, imageRegion4.GetSize()[1]); imageRegion3.SetSize(2, imageRegion4.GetSize()[2]); m_QballImage = QballImageType::New(); m_QballImage->SetSpacing( spacing3 ); m_QballImage->SetOrigin( origin3 ); m_QballImage->SetDirection( direction3 ); m_QballImage->SetRegions( imageRegion3 ); m_QballImage->Allocate(); Vector< PixelType, QBALL_ODFSIZE > nullVec1; nullVec1.Fill(0.0); m_QballImage->FillBuffer(nullVec1); m_CoefficientImage = CoefficientImageType::New(); m_CoefficientImage->SetSpacing( spacing3 ); m_CoefficientImage->SetOrigin( origin3 ); m_CoefficientImage->SetDirection( direction3 ); m_CoefficientImage->SetRegions( imageRegion3 ); m_CoefficientImage->Allocate(); Vector< PixelType, (ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder > nullVec2; nullVec2.Fill(0.0); m_CoefficientImage->FillBuffer(nullVec2); int x = imageRegion4.GetSize(0); int y = imageRegion4.GetSize(1); int z = imageRegion4.GetSize(2); int numCoeffs = imageRegion4.GetSize(3); for (int a=0; a coeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder,1); typename InputImageType::IndexType index; index.SetElement(0,a); index.SetElement(1,b); index.SetElement(2,c); typename CoefficientImageType::PixelType pix; for (int d=0; dGetPixel(index); coeffs[d][0] = pix[d]; } typename CoefficientImageType::IndexType index2; index2.SetElement(0,a); index2.SetElement(1,b); index2.SetElement(2,c); m_CoefficientImage->SetPixel(index2, pix); typename QballImageType::PixelType pix2; vnl_matrix odf = m_ShBasis*coeffs; for (int d=0; dSetPixel(index2,pix2); } } // generate spherical harmonic values of the desired order for each input direction template< class PixelType, int ShOrder > void ShCoefficientImageImporter< PixelType, ShOrder > ::CalcShBasis() { vnl_matrix_fixed sphCoords = GetSphericalOdfDirections(); int j, m; double mag, plm; for (int p=0; p(l,abs(m),cos(sphCoords(0,p))); mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; if (m<0) m_ShBasis(p,j) = sqrt(2.0)*mag*cos(-m*sphCoords(1,p)); else if (m==0) m_ShBasis(p,j) = mag; else m_ShBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(1,p)); break; case MRTRIX: plm = legendre_p(l,abs(m),-cos(sphCoords(0,p))); mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; if (m>0) m_ShBasis(p,j) = mag*cos(m*sphCoords(1,p)); else if (m==0) m_ShBasis(p,j) = mag; else m_ShBasis(p,j) = mag*sin(-m*sphCoords(1,p)); break; } j++; } } } // convert cartesian to spherical coordinates template< class PixelType, int ShOrder > vnl_matrix_fixed ShCoefficientImageImporter< PixelType, ShOrder > ::GetSphericalOdfDirections() { itk::OrientationDistributionFunction< PixelType, QBALL_ODFSIZE > odf; vnl_matrix_fixed* dir = odf.GetDirections(); vnl_matrix_fixed sphCoords; for (int i=0; iget_column(i).magnitude(); if( magget(2,i)/mag); // theta sphCoords(1,i) = atan2(dir->get(1,i), dir->get(0,i)); // phi } } return sphCoords; } } #endif // __itkShCoefficientImageImporter_cpp diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp index ec2d6b64e9..397dd300d9 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerOdf.cpp @@ -1,381 +1,275 @@ /*=================================================================== 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 "mitkTrackingHandlerOdf.h" #include #include #include #include +#include namespace mitk { TrackingHandlerOdf::TrackingHandlerOdf() - : m_GfaThreshold(0.1) - , m_SecondOrder(false) - , m_MinMaxNormalize(true) - , m_NumProbSamples(10) + : m_GfaThreshold(0.2) + , m_OdfThreshold(0.1) + , m_SharpenOdfs(false) + , m_NumProbSamples(1) { } TrackingHandlerOdf::~TrackingHandlerOdf() { } void TrackingHandlerOdf::InitForTracking() { MITK_INFO << "Initializing ODF tracker."; if (m_NeedsDataInit) { m_OdfHemisphereIndices.clear(); m_OdfReducedIndices.clear(); itk::OrientationDistributionFunction< float, QBALL_ODFSIZE > odf; vnl_vector_fixed ref; ref.fill(0); ref[0]=1; for (int i=0; i0) m_OdfHemisphereIndices.push_back(i); m_OdfFloatDirs.set_size(m_OdfHemisphereIndices.size(), 3); for (unsigned int i=0; i GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(m_OdfImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); m_GfaImage = gfaFilter->GetOutput(); } -// m_WorkingOdfImage = ItkOdfImageType::New(); -// m_WorkingOdfImage->SetSpacing( m_OdfImage->GetSpacing() ); -// m_WorkingOdfImage->SetOrigin( m_OdfImage->GetOrigin() ); -// m_WorkingOdfImage->SetDirection( m_OdfImage->GetDirection() ); -// m_WorkingOdfImage->SetRegions( m_OdfImage->GetLargestPossibleRegion() ); -// m_WorkingOdfImage->Allocate(); - -// MITK_INFO << "Rescaling ODFs."; -// ItkOdfImageType::SizeType image_size = m_OdfImage->GetLargestPossibleRegion().GetSize(); -//#ifdef WIN32 -//#pragma omp parallel for -//#else -//#pragma omp parallel for collapse(3) -//#endif -// for (unsigned int x=0; xGetPixel(idx); - -// for (int i=0; i=max) -// max = odf_values[i]; -// else if (odf_values[i]0.0) -// { -// odf_values -= min; -// odf_values /= max; -// } -// } -// else if (sum>0) -// odf_values /= sum; - -//#pragma omp critical -// m_WorkingOdfImage->SetPixel(idx, odf_values); -// } - m_NeedsDataInit = false; } + + std::cout << "TrackingHandlerOdf - GFA threshold: " << m_GfaThreshold << std::endl; + std::cout << "TrackingHandlerOdf - ODF threshold: " << m_OdfThreshold << std::endl; + if (m_SharpenOdfs) + std::cout << "TrackingHandlerOdf - Sharpening ODfs" << std::endl; } -vnl_vector< float > TrackingHandlerOdf::GetSecondOrderProbabilities(const itk::Point& itkP, vnl_vector< float >& angles, vnl_vector< float >& probs) +int TrackingHandlerOdf::SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles) { - vnl_vector< float > out_probs; - out_probs.set_size(m_OdfHemisphereIndices.size()); - out_probs.fill(0.0); - float out_probs_sum = 0; + boost::random::discrete_distribution dist(probs.begin(), probs.end()); + int sampled_idx = 0; + int max_sample_idx = -1; + float max_prob = 0; + int trials = 0; - int c = 0; - for (unsigned int j=0; j=m_AngularThreshold) - continue; - - vnl_vector_fixed d = m_OdfFloatDirs.get_row(j); - itk::Point pos; - pos[0] = itkP[0] + d[0]; - pos[1] = itkP[1] + d[1]; - pos[2] = itkP[2] + d[2]; - ItkOdfImageType::PixelType odf_values = GetImageValue(pos, m_OdfImage, m_Interpolate); - vnl_vector< float > new_angles = m_OdfFloatDirs*d; - - float probs_sum = 0; - vnl_vector< float > new_probs; new_probs.set_size(m_OdfHemisphereIndices.size()); - for (unsigned int i=0; i=m_AngularThreshold) - { - new_probs[i] = odf_values[m_OdfHemisphereIndices[i]]; - probs_sum += new_probs[i]; - } - else - new_probs[i] = 0; + boost::random::variate_generator> sampler(m_Rng, dist); + sampled_idx = sampler(); } - if (probs_sum>0.0001) - new_probs /= probs_sum; - - for (unsigned int i=0; imax_prob && probs[sampled_idx]>m_OdfThreshold && fabs(angles[sampled_idx])>=m_AngularThreshold) { - float p = new_probs[i] * probs[i]; - out_probs_sum += p; - out_probs[i] += p; + max_prob = probs[sampled_idx]; + max_sample_idx = sampled_idx; } - c += 1; + else if ( (probs[sampled_idx]<=m_OdfThreshold || fabs(angles[sampled_idx])0.0001) - out_probs /= out_probs_sum; - - return out_probs; -} - -bool TrackingHandlerOdf::MinMaxNormalize() const -{ - return m_MinMaxNormalize; -} - -void TrackingHandlerOdf::SetMinMaxNormalize(bool MinMaxNormalize) -{ - m_MinMaxNormalize = MinMaxNormalize; -} - -void TrackingHandlerOdf::SetSecondOrder(bool SecondOrder) -{ - m_SecondOrder = SecondOrder; + return max_sample_idx; } vnl_vector_fixed TrackingHandlerOdf::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); itk::Index<3> idx; m_OdfImage->TransformPhysicalPointToIndex(pos, idx); if ( !m_OdfImage->GetLargestPossibleRegion().IsInside(idx) ) return output_direction; + // check GFA threshold for termination float gfa = GetImageValue(pos, m_GfaImage, m_Interpolate); if (gfa last_dir; if (!olddirs.empty()) last_dir = olddirs.back(); if (!m_Interpolate && oldIndex==idx) return last_dir; ItkOdfImageType::PixelType odf_values = GetImageValue(pos, m_OdfImage, m_Interpolate); vnl_vector< float > probs; probs.set_size(m_OdfHemisphereIndices.size()); + vnl_vector< float > angles; angles.set_size(m_OdfHemisphereIndices.size()); angles.fill(1.0); + // Find ODF maximum and remove <0 values float max_odf_val = 0; - float odf_sum = 0; + float min_odf_val = 999; int max_idx_d = -1; int c = 0; for (int i : m_OdfHemisphereIndices) { if (odf_values[i]<0) odf_values[i] = 0; + if (odf_values[i]>max_odf_val) { max_odf_val = odf_values[i]; max_idx_d = c; } + if (odf_values[i]=0 && (olddirs.empty() || last_dir.magnitude()<=0.5)) // no previous direction, so return principal diffusion direction + if (m_SharpenOdfs) { - if (m_Mode==MODE::DETERMINISTIC) + // sharpen ODF + probs -= min_odf_val; + probs /= (max_odf_val-min_odf_val); + for (unsigned int i=0; i0) + { + probs /= odf_sum; + max_odf_val /= odf_sum; + } + } + + // no previous direction + if (max_odf_val>m_OdfThreshold && (olddirs.empty() || last_dir.magnitude()<=0.5)) + { + if (m_Mode==MODE::DETERMINISTIC) // return maximum peak { output_direction = m_OdfFloatDirs.get_row(max_idx_d); return output_direction * max_odf_val; } - else if (m_Mode==MODE::PROBABILISTIC) + else if (m_Mode==MODE::PROBABILISTIC) // sample from complete ODF { - probs /= odf_sum; - boost::random::discrete_distribution dist(probs.begin(), probs.end()); - int sampled_idx = 0; - int max_sample_idx = -1; - float max_prob = 0; - int trials = 0; - - for (int i=0; i> sampler(m_Rng, dist); - sampled_idx = sampler(); - } - if (probs[sampled_idx]>max_prob && probs[sampled_idx]>1.2/QBALL_ODFSIZE) - { - max_prob = probs[sampled_idx]; - max_sample_idx = sampled_idx; - } - else if (probs[sampled_idx]<=1.2/QBALL_ODFSIZE && trials<50) - i--; - } + int max_sample_idx = SampleOdf(probs, angles); if (max_sample_idx>=0) - output_direction = m_OdfFloatDirs.get_row(max_sample_idx); + output_direction = m_OdfFloatDirs.get_row(max_sample_idx) * probs[max_sample_idx]; return output_direction; } } - else if (max_idx_d<0) + else if (max_odf_val<=m_OdfThreshold) // return (0,0,0) { return output_direction; } + // correct previous direction if (m_FlipX) last_dir[0] *= -1; if (m_FlipY) last_dir[1] *= -1; if (m_FlipZ) last_dir[2] *= -1; - vnl_vector< float > angles = m_OdfFloatDirs*last_dir; + // calculate angles between previous direction and ODF directions + angles = m_OdfFloatDirs*last_dir; + float probs_sum = 0; float max_prob = 0; for (unsigned int i=0; imax_prob) + if (m_Mode==MODE::DETERMINISTIC && odf_val>max_prob && odf_val>m_OdfThreshold) { + // use maximum peak of the ODF weighted with the directional prior max_prob = odf_val; vnl_vector_fixed d = m_OdfFloatDirs.get_row(i); - if (angle<0) // make sure we don't walk backwards + if (angle<0) d *= -1; output_direction = odf_val*d; } else if (m_Mode==MODE::PROBABILISTIC) { + // update ODF probabilties with the ODF values pow(abs_angle, m_DirPriorPower) probs[i] = odf_val; probs_sum += probs[i]; } } + + // do probabilistic sampling if (m_Mode==MODE::PROBABILISTIC && probs_sum>0.0001) { - probs /= probs_sum; - - if (m_SecondOrder) - probs = GetSecondOrderProbabilities(pos, angles, probs); - - boost::random::discrete_distribution dist(probs.begin(), probs.end()); - - int sampled_idx = 0; - int max_sample_idx = -1; - float max_prob = 0; - int trials = 0; - - for (int i=0; i> sampler(m_Rng, dist); - sampled_idx = sampler(); - } - if (probs[sampled_idx]>max_prob && probs[sampled_idx]>1.2/QBALL_ODFSIZE) - { - max_prob = probs[sampled_idx]; - max_sample_idx = sampled_idx; - } - else if (probs[sampled_idx]<=1.2/QBALL_ODFSIZE && trials<50) - i--; - } +// probs /= probs_sum; + int max_sample_idx = SampleOdf(probs, angles); if (max_sample_idx>=0) { output_direction = m_OdfFloatDirs.get_row(max_sample_idx); - if (angles[max_sample_idx]<0) // make sure we don't walk backwards + if (angles[max_sample_idx]<0) output_direction *= -1; - output_direction *= max_prob; + output_direction *= probs[max_sample_idx]; } } + // check hard angular threshold float mag = output_direction.magnitude(); if (mag>=0.0001) { output_direction.normalize(); float a = dot_product(output_direction, last_dir); if (a #include #include namespace mitk { /** * \brief Enables streamline tracking on tensor images. Supports multi tensor tracking by adding multiple tensor images. */ class MITKFIBERTRACKING_EXPORT TrackingHandlerOdf : public TrackingDataHandler { public: TrackingHandlerOdf(); ~TrackingHandlerOdf(); typedef itk::DiffusionTensor3D TensorType; typedef itk::Image< itk::Vector< float, QBALL_ODFSIZE >, 3 > ItkOdfImageType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; void InitForTracking(); ///< calls InputDataValidForTracking() and creates feature images vnl_vector_fixed ProposeDirection(const itk::Point& pos, std::deque< vnl_vector_fixed >& olddirs, itk::Index<3>& oldIndex); ///< predicts next progression direction at the given position - + void SetSharpenOdfs(bool doSharpen) { m_SharpenOdfs=doSharpen; } + void SetOdfThreshold(float odfThreshold){ m_OdfThreshold = odfThreshold; } void SetGfaThreshold(float gfaThreshold){ m_GfaThreshold = gfaThreshold; } void SetOdfImage( ItkOdfImageType::Pointer img ){ m_OdfImage = img; DataModified(); } void SetGfaImage( ItkFloatImgType::Pointer img ){ m_GfaImage = img; DataModified(); } void SetMode( MODE m ){ m_Mode = m; } ItkUcharImgType::SpacingType GetSpacing(){ return m_OdfImage->GetSpacing(); } itk::Point GetOrigin(){ return m_OdfImage->GetOrigin(); } ItkUcharImgType::DirectionType GetDirection(){ return m_OdfImage->GetDirection(); } ItkUcharImgType::RegionType GetLargestPossibleRegion(){ return m_OdfImage->GetLargestPossibleRegion(); } int OdfPower() const; void SetNumProbSamples(int NumProbSamples); - void SetSecondOrder(bool SecondOrder); - - bool MinMaxNormalize() const; - void SetMinMaxNormalize(bool MinMaxNormalize); - protected: - vnl_vector< float > GetSecondOrderProbabilities(const itk::Point& itkP, vnl_vector< float >& angles, vnl_vector< float >& probs); + int SampleOdf(vnl_vector< float >& probs, vnl_vector< float >& angles); - float m_GfaThreshold; + float m_GfaThreshold; + float m_OdfThreshold; + bool m_SharpenOdfs; ItkFloatImgType::Pointer m_GfaImage; ///< GFA image used to determine streamline termination. ItkOdfImageType::Pointer m_OdfImage; ///< Input odf image. ItkOdfImageType::Pointer m_WorkingOdfImage; ///< Modified odf image. std::vector< int > m_OdfHemisphereIndices; vnl_matrix< float > m_OdfFloatDirs; - bool m_SecondOrder; - bool m_MinMaxNormalize; int m_NumProbSamples; std::vector< int > m_OdfReducedIndices; }; } #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp index 8ae378a920..f9ae39c147 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerPeaks.cpp @@ -1,279 +1,281 @@ /*=================================================================== 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) { std::srand(0); } 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 (magmitk::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/TrackingHandlers/mitkTrackingHandlerTensor.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp index 7960673741..3de948d392 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/TrackingHandlers/mitkTrackingHandlerTensor.cpp @@ -1,368 +1,372 @@ /*=================================================================== 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 "mitkTrackingHandlerTensor.h" namespace mitk { TrackingHandlerTensor::TrackingHandlerTensor() : m_InterpolateTensors(true) , m_FaThreshold(0.1) { } TrackingHandlerTensor::~TrackingHandlerTensor() { } void TrackingHandlerTensor::InitForTracking() { MITK_INFO << "Initializing tensor tracker."; if (m_NeedsDataInit) { m_NumberOfInputs = m_TensorImages.size(); for (int i=0; iSetSpacing( m_TensorImages.at(0)->GetSpacing() ); pdImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); pdImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); pdImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); pdImage->Allocate(); m_PdImage.push_back(pdImage); ItkDoubleImgType::Pointer emaxImage = ItkDoubleImgType::New(); emaxImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); emaxImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); emaxImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); emaxImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); emaxImage->Allocate(); emaxImage->FillBuffer(1.0); m_EmaxImage.push_back(emaxImage); } bool useUserFaImage = true; if (m_FaImage.IsNull()) { m_FaImage = ItkFloatImgType::New(); m_FaImage->SetSpacing( m_TensorImages.at(0)->GetSpacing() ); m_FaImage->SetOrigin( m_TensorImages.at(0)->GetOrigin() ); m_FaImage->SetDirection( m_TensorImages.at(0)->GetDirection() ); m_FaImage->SetRegions( m_TensorImages.at(0)->GetLargestPossibleRegion() ); m_FaImage->Allocate(); m_FaImage->FillBuffer(0.0); useUserFaImage = false; } typedef itk::DiffusionTensor3D TensorType; TensorType::EigenValuesArrayType eigenvalues; TensorType::EigenVectorsMatrixType eigenvectors; #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif for (int x=0; xGetLargestPossibleRegion().GetSize()[0]; x++) for (int y=0; yGetLargestPossibleRegion().GetSize()[1]; y++) for (int z=0; zGetLargestPossibleRegion().GetSize()[2]; z++) { ItkTensorImageType::IndexType index; index[0] = x; index[1] = y; index[2] = z; for (int i=0; iGetPixel(index); tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); } vnl_vector_fixed dir; dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude()>mitk::eps) dir.normalize(); else dir.fill(0.0); #pragma omp critical { m_PdImage.at(i)->SetPixel(index, dir); if (!useUserFaImage) m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)+tensor.GetFractionalAnisotropy()); m_EmaxImage.at(i)->SetPixel(index, 2/eigenvalues[2]); } } if (!useUserFaImage) { #pragma omp critical m_FaImage->SetPixel(index, m_FaImage->GetPixel(index)/m_NumberOfInputs); } } m_NeedsDataInit = false; } if (m_F+m_G>1.0) { float temp = m_F+m_G; m_F /= temp; m_G /= temp; } + + std::cout << "TrackingHandlerTensor - FA threshold: " << m_FaThreshold << std::endl; + std::cout << "TrackingHandlerTensor - f: " << m_F << std::endl; + std::cout << "TrackingHandlerTensor - g: " << m_G << std::endl; } vnl_vector_fixed TrackingHandlerTensor::GetMatchingDirection(itk::Index<3> idx, vnl_vector_fixed& oldDir, int& image_num) { vnl_vector_fixed out_dir; out_dir.fill(0); float angle = 0; float mag = oldDir.magnitude(); if (magGetPixel(idx); if (out_dir.magnitude()>0.5) { image_num = i; oldDir[0] = out_dir[0]; oldDir[1] = out_dir[1]; oldDir[2] = out_dir[2]; break; } } } else { for (unsigned int i=0; i dir = m_PdImage.at(i)->GetPixel(idx); float a = dot_product(dir, oldDir); if (fabs(a)>angle) { image_num = i; angle = fabs(a); if (a<0) out_dir = -dir; else out_dir = dir; out_dir *= angle; // shrink contribution of direction if is less parallel to previous direction } } } return out_dir; } vnl_vector_fixed TrackingHandlerTensor::GetDirection(itk::Point itkP, vnl_vector_fixed oldDir, TensorType& tensor) { // transform physical point to index coordinates itk::Index<3> idx; itk::ContinuousIndex< float, 3> cIdx; m_FaImage->TransformPhysicalPointToIndex(itkP, idx); m_FaImage->TransformPhysicalPointToContinuousIndex(itkP, cIdx); vnl_vector_fixed dir; dir.fill(0.0); if ( !m_FaImage->GetLargestPossibleRegion().IsInside(idx) ) return dir; int image_num = -1; if (!m_Interpolate) { dir = GetMatchingDirection(idx, oldDir, image_num); if (image_num>=0) tensor = m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx); } else { float frac_x = cIdx[0] - idx[0]; float frac_y = cIdx[1] - idx[1]; float frac_z = cIdx[2] - idx[2]; if (frac_x<0) { idx[0] -= 1; frac_x += 1; } if (frac_y<0) { idx[1] -= 1; frac_y += 1; } if (frac_z<0) { idx[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 (idx[0] >= 0 && idx[0] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(0) - 1) && idx[1] >= 0 && idx[1] < static_cast(m_FaImage->GetLargestPossibleRegion().GetSize(1) - 1) && idx[2] >= 0 && idx[2] < static_cast(m_FaImage->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(idx, oldDir, image_num) * interpWeights[0]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(idx) * m_EmaxImage[image_num]->GetPixel(idx) * interpWeights[0]; itk::Index<3> tmpIdx = idx; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[1]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = idx; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[2]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = idx; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[3]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[4]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = idx; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[5]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = idx; tmpIdx[2]++; tmpIdx[0]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[6]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = idx; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; dir += GetMatchingDirection(tmpIdx, oldDir, image_num) * interpWeights[7]; if (image_num>=0) tensor += m_TensorImages[image_num]->GetPixel(tmpIdx) * m_EmaxImage[image_num]->GetPixel(tmpIdx) * interpWeights[7]; } } return dir; } vnl_vector_fixed TrackingHandlerTensor::GetLargestEigenvector(TensorType& tensor) { vnl_vector_fixed dir; TensorType::EigenValuesArrayType eigenvalues; TensorType::EigenVectorsMatrixType eigenvectors; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); if (dir.magnitude() TrackingHandlerTensor::ProposeDirection(const itk::Point& pos, std::deque >& olddirs, itk::Index<3>& oldIndex) { vnl_vector_fixed output_direction; output_direction.fill(0); TensorType tensor; tensor.Fill(0); try { itk::Index<3> index; m_TensorImages.at(0)->TransformPhysicalPointToIndex(pos, index); float fa = GetImageValue(pos, m_FaImage, m_Interpolate); if (fa oldDir = olddirs.back(); if (m_FlipX) oldDir[0] *= -1; if (m_FlipY) oldDir[1] *= -1; if (m_FlipZ) oldDir[2] *= -1; float old_mag = oldDir.magnitude(); if (!m_Interpolate && oldIndex==index) return oldDir; output_direction = GetDirection(pos, oldDir, tensor); float mag = output_direction.magnitude(); if (mag>=mitk::eps) { output_direction.normalize(); if (old_mag>0.5 && m_G>mitk::eps) // TEND tracking { output_direction[0] = m_F*output_direction[0] + (1-m_F)*( (1-m_G)*oldDir[0] + m_G*(tensor[0]*oldDir[0] + tensor[1]*oldDir[1] + tensor[2]*oldDir[2])); output_direction[1] = m_F*output_direction[1] + (1-m_F)*( (1-m_G)*oldDir[1] + m_G*(tensor[1]*oldDir[0] + tensor[3]*oldDir[1] + tensor[4]*oldDir[2])); output_direction[2] = m_F*output_direction[2] + (1-m_F)*( (1-m_G)*oldDir[2] + m_G*(tensor[2]*oldDir[0] + tensor[4]*oldDir[1] + tensor[5]*oldDir[2])); 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); } catch(...) { } if (m_FlipX) output_direction[0] *= -1; if (m_FlipY) output_direction[1] *= -1; if (m_FlipZ) output_direction[2] *= -1; return output_direction; } } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp index 944683e738..6763b1be91 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkStreamlineTrackingFilter.cpp @@ -1,1010 +1,1010 @@ /*=================================================================== 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 #define _USE_MATH_DEFINES #include namespace itk { StreamlineTrackingFilter ::StreamlineTrackingFilter() : m_PauseTracking(false) , m_AbortTracking(false) , m_FiberPolyData(nullptr) , m_Points(nullptr) , m_Cells(nullptr) , m_AngularThreshold(-1) , m_StepSize(0) , m_MaxLength(10000) , m_MinTractLength(20.0) , m_MaxTractLength(400.0) , m_SeedsPerVoxel(1) , m_RandomSampling(false) , m_SamplingDistance(-1) , m_DeflectionMod(1.0) , m_OnlyForwardSamples(true) , m_UseStopVotes(true) , m_NumberOfSamples(30) , m_NumPreviousDirections(1) , m_StoppingRegions(nullptr) , m_SeedImage(nullptr) , m_MaskImage(nullptr) , m_AposterioriCurvCheck(false) , m_AvoidStop(true) , m_DemoMode(false) , m_SeedOnlyGm(false) , m_ControlGmEndings(false) , m_WmLabel(3) // mrtrix 5ttseg labels , m_GmLabel(1) // mrtrix 5ttseg labels , m_StepSizeVox(-1) , m_SamplingDistanceVox(-1) , m_AngularThresholdDeg(-1) , m_MaxNumTracts(-1) , m_Random(true) , m_Verbose(true) , m_UseOutputProbabilityMap(false) { this->SetNumberOfRequiredInputs(0); } void StreamlineTrackingFilter::BeforeTracking() { 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; } float StreamlineTrackingFilter::GetRandDouble(float min, float max) { return (float)(rand()%((int)(10000*(max-min))) + 10000*min)/10000; } std::vector< vnl_vector_fixed > StreamlineTrackingFilter::CreateDirections(int NPoints) { std::vector< vnl_vector_fixed > pointshell; if (NPoints<2) return pointshell; std::vector< float > theta; theta.resize(NPoints); std::vector< float > phi; phi.resize(NPoints); float C = sqrt(4*M_PI); phi[0] = 0.0; phi[NPoints-1] = 0.0; for(int i=0; i0 && i d; d[0] = cos(theta[i]) * cos(phi[i]); d[1] = cos(theta[i]) * sin(phi[i]); d[2] = sin(theta[i]); pointshell.push_back(d); } return pointshell; } vnl_vector_fixed StreamlineTrackingFilter::GetNewDirection(itk::Point &pos, std::deque >& olddirs, itk::Index<3> &oldIndex) { if (m_DemoMode) { m_SamplingPointset->Clear(); m_AlternativePointset->Clear(); m_StopVotePointset->Clear(); } vnl_vector_fixed direction; direction.fill(0); ItkUcharImgType::IndexType idx; m_MaskImage->TransformPhysicalPointToIndex(pos, idx); if (IsValidPosition(pos)) { if (m_StoppingRegions->GetPixel(idx)>0) return direction; direction = m_TrackingHandler->ProposeDirection(pos, olddirs, oldIndex); // get direction proposal at current streamline position } else return direction; vnl_vector_fixed olddir = olddirs.back(); std::vector< vnl_vector_fixed > probeVecs = CreateDirections(m_NumberOfSamples); itk::Point sample_pos; int alternatives = 1; int stop_votes = 0; int possible_stop_votes = 0; for (int i=0; i d; bool is_stop_voter = false; if (m_RandomSampling) { d[0] = GetRandDouble(); d[1] = GetRandDouble(); d[2] = GetRandDouble(); d.normalize(); d *= GetRandDouble(0,m_SamplingDistance); } else { d = probeVecs.at(i); float dot = dot_product(d, olddir); if (m_UseStopVotes && dot>0.7) { is_stop_voter = true; possible_stop_votes++; } else if (m_OnlyForwardSamples && dot<0) continue; d *= m_SamplingDistance; } sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (IsValidPosition(sample_pos)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) { direction += tempDir; if(m_DemoMode) m_SamplingPointset->InsertPoint(i, sample_pos); } else if (m_AvoidStop && olddir.magnitude()>0.5) // out of white matter { if (is_stop_voter) stop_votes++; if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); float dot = dot_product(d, olddir); if (dot >= 0.0) // in front of plane defined by pos and olddir d = -d + 2*dot*olddir; // reflect else d = -d; // invert // look a bit further into the other direction sample_pos[0] = pos[0] + d[0]; sample_pos[1] = pos[1] + d[1]; sample_pos[2] = pos[2] + d[2]; alternatives++; vnl_vector_fixed tempDir; tempDir.fill(0.0); if (IsValidPosition(sample_pos)) tempDir = m_TrackingHandler->ProposeDirection(sample_pos, olddirs, oldIndex); // sample neighborhood if (tempDir.magnitude()>mitk::eps) // are we back in the white matter? { direction += d * m_DeflectionMod; // go into the direction of the white matter direction += tempDir; // go into the direction of the white matter direction at this location if(m_DemoMode) m_AlternativePointset->InsertPoint(alternatives, sample_pos); } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); } } else { if (m_DemoMode) m_StopVotePointset->InsertPoint(i, sample_pos); if (is_stop_voter) stop_votes++; } } if (direction.magnitude()>0.001 && (possible_stop_votes==0 || (float)stop_votes/possible_stop_votes<0.5) ) direction.normalize(); else direction.fill(0); return direction; } float StreamlineTrackingFilter::FollowStreamline(itk::Point pos, vnl_vector_fixed dir, FiberType* fib, float tractLength, bool front) { vnl_vector_fixed zero_dir; zero_dir.fill(0.0); std::deque< vnl_vector_fixed > last_dirs; for (unsigned int i=0; iTransformPhysicalPointToIndex(pos, oldIndex); // get new position CalculateNewPosition(pos, dir); // is new position inside of image and mask if (m_AbortTracking) // if not end streamline { return tractLength; } else // if yes, add new point to streamline { tractLength += m_StepSize; if (front) fib->push_front(pos); else fib->push_back(pos); if (m_AposterioriCurvCheck) { int curv = CheckCurvature(fib, front); if (curv>0) { tractLength -= m_StepSize*curv; while (curv>0) { if (front) fib->pop_front(); else fib->pop_back(); curv--; } return tractLength; } } if (tractLength>m_MaxTractLength) return tractLength; } 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==fib->size()-1) meanV += v; c--; } } meanV.normalize(); for (int c=0; c1.0) angle = 1.0; if (angle<-1.0) angle = -1.0; dev += acos(angle)*180/M_PI; } if (vectors.size()>0) dev /= vectors.size(); if (dev<30) return 0; else return vectors.size(); } void StreamlineTrackingFilter::GetSeedPointsFromSeedImage() { - MITK_INFO << "Calculating seed points."; + 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] + GetRandDouble(-0.5, 0.5); start[1] = index[1] + GetRandDouble(-0.5, 0.5); start[2] = index[2] + GetRandDouble(-0.5, 0.5); itk::Point worldPos; m_SeedImage->TransformContinuousIndexToPhysicalPoint(start, worldPos); m_SeedPoints.push_back(worldPos); } } } ++sit; } } else { for (auto s : m_GmStubs) m_SeedPoints.push_back(s[1]); } } void StreamlineTrackingFilter::GenerateData() { this->BeforeTracking(); if (m_Random) { std::srand(std::time(0)); std::random_shuffle(m_SeedPoints.begin(), m_SeedPoints.end()); } bool stop = false; unsigned int current_tracts = 0; int num_seeds = m_SeedPoints.size(); itk::Index<3> zeroIndex; zeroIndex.Fill(0); int progress = 0; int i = 0; int print_interval = num_seeds/100; if (print_interval<100) m_Verbose=false; #pragma omp parallel while (i=num_seeds || stop) continue; else if (m_Verbose && i%print_interval==0) #pragma omp critical { progress += print_interval; std::cout << " \r"; if (m_MaxNumTracts>0) std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << "/" << m_MaxNumTracts << '\r'; else std::cout << "Tried: " << progress << "/" << num_seeds << " | Accepted: " << current_tracts << '\r'; cout.flush(); } 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_UseOutputProbabilityMap) m_Tractogram.push_back(fib); else FiberToProbmap(&fib); current_tracts++; } if (m_MaxNumTracts > 0 && current_tracts>=static_cast(m_MaxNumTracts)) { if (!stop) { std::cout << " \r"; MITK_INFO << "Reconstructed maximum number of tracts (" << current_tracts << "). Stopping tractography."; } stop = true; } } } } this->AfterTracking(); } void StreamlineTrackingFilter::CheckFiberForGmEnding(FiberType* fib) { if (m_TissueImage.IsNull()) return; // first check if the current fibe rendpoint is located inside of the white matter // if not, remove last fiber point and repeat bool in_wm = false; while (!in_wm && fib->size()>2) { ItkUcharImgType::IndexType idx; m_TissueImage->TransformPhysicalPointToIndex(fib->back(), idx); if (m_TissueImage->GetPixel(idx)==m_WmLabel) in_wm = true; else fib->pop_back(); } if (fib->size()<3 || !in_wm) { fib->clear(); return; } // get fiber direction at end point vnl_vector_fixed< float, 3 > d1; d1[0] = fib->back()[0] - fib->at(fib->size()-2)[0]; d1[1] = fib->back()[1] - fib->at(fib->size()-2)[1]; d1[2] = fib->back()[2] - fib->at(fib->size()-2)[2]; d1.normalize(); // find closest gray matter voxel ItkUcharImgType::IndexType s_idx; m_TissueImage->TransformPhysicalPointToIndex(fib->back(), s_idx); itk::Point gm_endp; float max = -1; for (int x : {-1,0,1}) for (int y : {-1,0,1}) for (int z : {-1,0,1}) { if (x==y && y==z) continue; ItkUcharImgType::IndexType e_idx; e_idx[0] = s_idx[0] + x; e_idx[1] = s_idx[1] + y; e_idx[2] = s_idx[2] + z; if ( !m_TissueImage->GetLargestPossibleRegion().IsInside(e_idx) || m_TissueImage->GetPixel(e_idx)!=m_GmLabel ) continue; itk::ContinuousIndex end; m_TissueImage->TransformIndexToPhysicalPoint(e_idx, end); vnl_vector_fixed< float, 3 > d2; d2[0] = end[0] - fib->back()[0]; d2[1] = end[1] - fib->back()[1]; d2[2] = end[2] - fib->back()[2]; d2.normalize(); float a = dot_product(d1,d2); if (a>max) { max = a; gm_endp = end; } } if (max>=0) fib->push_back(gm_endp); else // no gray matter enpoint found -> delete fiber fib->clear(); } void StreamlineTrackingFilter::FiberToProbmap(FiberType* fib) { ItkDoubleImgType::IndexType last_idx; last_idx.Fill(0); for (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 (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/cmdapps/Tractography/StreamlineTractography.cpp b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp index c843d43421..10e7331fc2 100755 --- a/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp +++ b/Modules/DiffusionImaging/FiberTracking/cmdapps/Tractography/StreamlineTractography.cpp @@ -1,466 +1,463 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #define _USE_MATH_DEFINES #include using namespace std; const int numOdfSamples = 200; typedef itk::Image< itk::Vector< float, numOdfSamples > , 3 > SampledShImageType; /*! \brief */ int main(int argc, char* argv[]) { mitkCommandLineParser parser; parser.setTitle("Streamline Tractography"); parser.setCategory("Fiber Tracking and Processing Methods"); parser.setDescription("Perform streamline tractography"); parser.setContributor("MIC"); // parameters fo all methods parser.setArgumentPrefix("--", "-"); parser.addArgument("input", "i", mitkCommandLineParser::StringList, "Input:", "input image (multiple possible for 'DetTensor' algorithm)", us::Any(), false); parser.addArgument("algorithm", "a", mitkCommandLineParser::String, "Algorithm:", "which algorithm to use (Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF)", us::Any(), false); parser.addArgument("out", "o", mitkCommandLineParser::OutputDirectory, "Output:", "output fiberbundle/probability map", us::Any(), false); parser.addArgument("stop_mask", "", mitkCommandLineParser::String, "Stop image:", "streamlines entering the binary mask will stop immediately", us::Any()); parser.addArgument("tracking_mask", "", mitkCommandLineParser::String, "Mask image:", "restrict tractography with a binary mask image", us::Any()); parser.addArgument("seed_mask", "", mitkCommandLineParser::String, "Seed image:", "binary mask image defining seed voxels", us::Any()); parser.addArgument("tissue_image", "", mitkCommandLineParser::String, "Tissue type image:", "image with tissue type labels (WM=3, GM=1)", us::Any()); + parser.addArgument("sharpen_odfs", "", mitkCommandLineParser::Bool, "SHarpen ODFs:", "if you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). this is not necessary for CSD fODFs, since they are narurally much sharper."); parser.addArgument("cutoff", "", mitkCommandLineParser::Float, "Cutoff:", "set the FA, GFA or Peak amplitude cutoff for terminating tracks", 0.1); + parser.addArgument("odf_cutoff", "", mitkCommandLineParser::Float, "ODF Cutoff:", "additional threshold on the ODF magnitude. this is useful in case of CSD fODF tractography.", 0.1); parser.addArgument("step_size", "", mitkCommandLineParser::Float, "Step size:", "step size (in voxels)", 0.5); parser.addArgument("angular_threshold", "", mitkCommandLineParser::Float, "Angular threshold:", "angular threshold between two successive steps, (default: 90° * step_size)"); parser.addArgument("seeds", "", mitkCommandLineParser::Int, "Seeds per voxel:", "number of seed points per voxel", 1); parser.addArgument("seed_gm", "", mitkCommandLineParser::Bool, "Seed only GM:", "Seed only in gray matter (requires tissue type image --tissue_image)"); parser.addArgument("control_gm_endings", "", mitkCommandLineParser::Bool, "Control GM endings:", "Seed perpendicular to gray matter and enforce endings inside of the gray matter (requires tissue type image --tissue_image)"); parser.addArgument("max_tracts", "", mitkCommandLineParser::Int, "Max. number of tracts:", "tractography is stopped if the reconstructed number of tracts is exceeded.", -1); parser.addArgument("num_samples", "", mitkCommandLineParser::Int, "Num. neighborhood samples:", "number of neighborhood samples that are use to determine the next progression direction", 0); parser.addArgument("sampling_distance", "", mitkCommandLineParser::Float, "Sampling distance:", "distance of neighborhood sampling points (in voxels)", 0.25); parser.addArgument("use_stop_votes", "", mitkCommandLineParser::Bool, "Use stop votes:", "use stop votes"); parser.addArgument("use_only_forward_samples", "", mitkCommandLineParser::Bool, "Use only forward samples:", "use only forward samples"); parser.addArgument("output_prob_map", "", mitkCommandLineParser::Bool, "Output probability map:", "output probability map instead of tractogram"); - parser.addArgument("prob_samples", "", mitkCommandLineParser::Int, "Num. Samples per Step:", "Only relevant for probabilistic tensor and ODF tractography. At each integration step the ODF/tensor is sampled n times and the most probable direction is used. The higher the number, the closer the behaviour will be to the deterministic ODF tractography.", 10); parser.addArgument("no_interpolation", "", mitkCommandLineParser::Bool, "Don't interpolate:", "don't interpolate image values"); parser.addArgument("flip_x", "", mitkCommandLineParser::Bool, "Flip X:", "multiply x-coordinate of direction proposal by -1"); parser.addArgument("flip_y", "", mitkCommandLineParser::Bool, "Flip Y:", "multiply y-coordinate of direction proposal by -1"); parser.addArgument("flip_z", "", mitkCommandLineParser::Bool, "Flip Z:", "multiply z-coordinate of direction proposal by -1"); //parser.addArgument("apply_image_rotation", "", mitkCommandLineParser::Bool, "Apply image rotation:", "applies image rotation to image peaks (only for 'Peaks' algorithm). This is necessary in some cases, e.g. if the peaks were obtained with MRtrix."); parser.addArgument("compress", "", mitkCommandLineParser::Float, "Compress:", "Compress output fibers using the given error threshold (in mm)"); parser.addArgument("additional_images", "", mitkCommandLineParser::StringList, "Additional images:", "specify a list of float images that hold additional information (FA, GFA, additional Features)", us::Any()); - // parameters for ODF based tractography - parser.addArgument("no_odf_normalization", "", mitkCommandLineParser::Bool, "Don't min-max normalize ODFs:", "No min-max normalization of ODFs"); - // parameters for random forest based tractography parser.addArgument("forest", "", mitkCommandLineParser::String, "Forest:", "input random forest (HDF5 file)", us::Any()); parser.addArgument("use_sh_features", "", mitkCommandLineParser::Bool, "Use SH features:", "use SH features"); // parameters for tensor tractography parser.addArgument("tend_f", "", mitkCommandLineParser::Float, "Weight f", "Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).", 1.0); parser.addArgument("tend_g", "", mitkCommandLineParser::Float, "Weight g", "Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)", 0.0); map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size()==0) return EXIT_FAILURE; mitkCommandLineParser::StringContainerType input_files = us::any_cast(parsedArgs["input"]); string outFile = us::any_cast(parsedArgs["out"]); string algorithm = us::any_cast(parsedArgs["algorithm"]); + bool sharpen_odfs = false; + if (parsedArgs.count("sharpen_odfs")) + sharpen_odfs = us::any_cast(parsedArgs["sharpen_odfs"]); + bool interpolate = true; if (parsedArgs.count("no_interpolation")) interpolate = !us::any_cast(parsedArgs["no_interpolation"]); bool use_sh_features = false; if (parsedArgs.count("use_sh_features")) use_sh_features = us::any_cast(parsedArgs["use_sh_features"]); bool seed_gm = false; if (parsedArgs.count("seed_gm")) seed_gm = us::any_cast(parsedArgs["seed_gm"]); bool control_gm_endings = false; if (parsedArgs.count("control_gm_endings")) control_gm_endings = us::any_cast(parsedArgs["control_gm_endings"]); bool use_stop_votes = false; if (parsedArgs.count("use_stop_votes")) use_stop_votes = us::any_cast(parsedArgs["use_stop_votes"]); bool use_only_forward_samples = false; if (parsedArgs.count("use_only_forward_samples")) use_only_forward_samples = us::any_cast(parsedArgs["use_only_forward_samples"]); bool output_prob_map = false; if (parsedArgs.count("output_prob_map")) output_prob_map = us::any_cast(parsedArgs["output_prob_map"]); bool flip_x = false; if (parsedArgs.count("flip_x")) flip_x = us::any_cast(parsedArgs["flip_x"]); bool flip_y = false; if (parsedArgs.count("flip_y")) flip_y = us::any_cast(parsedArgs["flip_y"]); bool flip_z = false; if (parsedArgs.count("flip_z")) flip_z = us::any_cast(parsedArgs["flip_z"]); bool apply_image_rotation = false; if (parsedArgs.count("apply_image_rotation")) apply_image_rotation = us::any_cast(parsedArgs["apply_image_rotation"]); float compress = -1; if (parsedArgs.count("compress")) compress = us::any_cast(parsedArgs["compress"]); string forestFile; if (parsedArgs.count("forest")) forestFile = us::any_cast(parsedArgs["forest"]); string maskFile = ""; if (parsedArgs.count("tracking_mask")) maskFile = us::any_cast(parsedArgs["tracking_mask"]); string seedFile = ""; if (parsedArgs.count("seed_mask")) seedFile = us::any_cast(parsedArgs["seed_mask"]); string stopFile = ""; if (parsedArgs.count("stop_mask")) stopFile = us::any_cast(parsedArgs["stop_mask"]); string tissueFile = ""; if (parsedArgs.count("tissue_image")) tissueFile = us::any_cast(parsedArgs["tissue_image"]); float cutoff = 0.1; if (parsedArgs.count("cutoff")) cutoff = us::any_cast(parsedArgs["cutoff"]); + float odf_cutoff = 0.1; + if (parsedArgs.count("odf_cutoff")) + odf_cutoff = us::any_cast(parsedArgs["odf_cutoff"]); + float stepsize = -1; if (parsedArgs.count("step_size")) stepsize = us::any_cast(parsedArgs["step_size"]); float sampling_distance = -1; if (parsedArgs.count("sampling_distance")) sampling_distance = us::any_cast(parsedArgs["sampling_distance"]); int num_samples = 0; if (parsedArgs.count("num_samples")) num_samples = us::any_cast(parsedArgs["num_samples"]); int seeds = 1; if (parsedArgs.count("seeds")) seeds = us::any_cast(parsedArgs["seeds"]); float tend_f = 1; if (parsedArgs.count("tend_f")) tend_f = us::any_cast(parsedArgs["tend_f"]); float tend_g = 0; if (parsedArgs.count("tend_g")) tend_g = us::any_cast(parsedArgs["tend_g"]); float angular_threshold = -1; if (parsedArgs.count("angular_threshold")) angular_threshold = us::any_cast(parsedArgs["angular_threshold"]); - bool normalize_odfs = true; - if (parsedArgs.count("no_odf_normalization")) - normalize_odfs = !us::any_cast(parsedArgs["no_odf_normalization"]); - unsigned int max_tracts = -1; if (parsedArgs.count("max_tracts")) max_tracts = us::any_cast(parsedArgs["max_tracts"]); - unsigned int prob_samples = 10; - if (parsedArgs.count("prob_samples")) - prob_samples = us::any_cast(parsedArgs["prob_samples"]); - if (prob_samples<=0) - prob_samples = 1; - // LOAD DATASETS mitkCommandLineParser::StringContainerType addFiles; if (parsedArgs.count("additional_images")) addFiles = us::any_cast(parsedArgs["additional_images"]); typedef itk::Image ItkUcharImgType; MITK_INFO << "loading input"; std::vector< mitk::Image::Pointer > input_images; for (unsigned int i=0; i(mitk::IOUtil::LoadImage(maskFile).GetPointer()); mask = ItkUcharImgType::New(); mitk::CastToItkImage(img, mask); } ItkUcharImgType::Pointer seed; if (!seedFile.empty()) { MITK_INFO << "loading seed image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(seedFile).GetPointer()); seed = ItkUcharImgType::New(); mitk::CastToItkImage(img, seed); } ItkUcharImgType::Pointer stop; if (!stopFile.empty()) { MITK_INFO << "loading stop image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(stopFile).GetPointer()); stop = ItkUcharImgType::New(); mitk::CastToItkImage(img, stop); } ItkUcharImgType::Pointer tissue; if (!tissueFile.empty()) { MITK_INFO << "loading tissue image"; mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(tissueFile).GetPointer()); tissue = ItkUcharImgType::New(); mitk::CastToItkImage(img, tissue); } MITK_INFO << "loading additional images"; typedef itk::Image ItkFloatImgType; std::vector< std::vector< ItkFloatImgType::Pointer > > addImages; addImages.push_back(std::vector< ItkFloatImgType::Pointer >()); for (auto file : addFiles) { mitk::Image::Pointer img = dynamic_cast(mitk::IOUtil::LoadImage(file).GetPointer()); ItkFloatImgType::Pointer itkimg = ItkFloatImgType::New(); mitk::CastToItkImage(img, itkimg); addImages.at(0).push_back(itkimg); } // ////////////////////////////////////////////////////////////////// // omp_set_num_threads(1); if (algorithm == "ProbTensor") { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkTensorImg = caster->GetOutput(); typedef itk::TensorImageToQBallImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkTensorImg ); filter->Update(); mitk::Image::Pointer image = mitk::Image::New(); FilterType::OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); input_images.clear(); input_images.push_back(image); algorithm = "ProbODF"; + + sharpen_odfs = true; + odf_cutoff = 0; } typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); mitk::TrackingDataHandler* handler; if (algorithm == "DetRF" || algorithm == "ProbRF") { if (use_sh_features) { handler = new mitk::TrackingHandlerRandomForest<6,28>(); dynamic_cast*>(handler)->LoadForest(forestFile); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } else { handler = new mitk::TrackingHandlerRandomForest<6,100>(); dynamic_cast*>(handler)->LoadForest(forestFile); dynamic_cast*>(handler)->AddDwi(input_images.at(0)); dynamic_cast*>(handler)->SetAdditionalFeatureImages(addImages); } if (algorithm == "ProbRF") handler->SetMode(mitk::TrackingDataHandler::MODE::PROBABILISTIC); } else if (algorithm == "Peaks") { handler = new mitk::TrackingHandlerPeaks(); typedef mitk::ImageToItk< mitk::TrackingHandlerPeaks::PeakImgType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_images.at(0)); caster->Update(); mitk::TrackingHandlerPeaks::PeakImgType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetPeakImage(itkImg); dynamic_cast(handler)->SetApplyDirectionMatrix(apply_image_rotation); dynamic_cast(handler)->SetPeakThreshold(cutoff); } else if (algorithm == "DetTensor") { handler = new mitk::TrackingHandlerTensor(); for (auto input_image : input_images) { typedef mitk::ImageToItk< mitk::TrackingHandlerTensor::ItkTensorImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_image); caster->Update(); mitk::TrackingHandlerTensor::ItkTensorImageType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->AddTensorImage(itkImg); } dynamic_cast(handler)->SetFaThreshold(cutoff); dynamic_cast(handler)->SetF(tend_f); dynamic_cast(handler)->SetG(tend_g); if (addImages.at(0).size()>0) dynamic_cast(handler)->SetFaImage(addImages.at(0).at(0)); } else if (algorithm == "DetODF" || algorithm == "ProbODF") { handler = new mitk::TrackingHandlerOdf(); for (auto input_image : input_images) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(input_image); caster->Update(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = caster->GetOutput(); dynamic_cast(handler)->SetOdfImage(itkImg); } dynamic_cast(handler)->SetGfaThreshold(cutoff); + dynamic_cast(handler)->SetOdfThreshold(odf_cutoff); + dynamic_cast(handler)->SetSharpenOdfs(sharpen_odfs); + if (algorithm == "ProbODF") - { dynamic_cast(handler)->SetMode(mitk::TrackingHandlerOdf::MODE::PROBABILISTIC); - dynamic_cast(handler)->SetNumProbSamples(prob_samples); - } if (addImages.at(0).size()>0) dynamic_cast(handler)->SetGfaImage(addImages.at(0).at(0)); - - dynamic_cast(handler)->SetMinMaxNormalize(normalize_odfs); } else { MITK_INFO << "Unknown tractography algorithm (" + algorithm+"). Known types are Peaks, DetTensor, ProbTensor, DetODF, ProbODF, DetRF, ProbRF."; return EXIT_FAILURE; } handler->SetInterpolate(interpolate); handler->SetFlipX(flip_x); handler->SetFlipY(flip_y); handler->SetFlipZ(flip_z); MITK_INFO << "Tractography algorithm: " << algorithm; tracker->SetNumberOfSamples(num_samples); tracker->SetAngularThreshold(angular_threshold); tracker->SetMaskImage(mask); tracker->SetSeedImage(seed); tracker->SetStoppingRegions(stop); tracker->SetSeedsPerVoxel(seeds); tracker->SetStepSize(stepsize); tracker->SetSamplingDistance(sampling_distance); tracker->SetUseStopVotes(use_stop_votes); tracker->SetOnlyForwardSamples(use_only_forward_samples); tracker->SetAposterioriCurvCheck(false); tracker->SetTissueImage(tissue); tracker->SetSeedOnlyGm(seed_gm); tracker->SetControlGmEndings(control_gm_endings); tracker->SetMaxNumTracts(max_tracts); tracker->SetTrackingHandler(handler); tracker->SetUseOutputProbabilityMap(output_prob_map); tracker->Update(); std::string ext = itksys::SystemTools::GetFilenameExtension(outFile); if (!output_prob_map) { vtkSmartPointer< vtkPolyData > poly = tracker->GetFiberPolyData(); mitk::FiberBundle::Pointer outFib = mitk::FiberBundle::New(poly); if (compress > 0) outFib->Compress(compress); if (ext != ".fib" && ext != ".trk" && ext != ".tck") outFile += ".fib"; mitk::IOUtil::Save(outFib, outFile); } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (ext != ".nii" && ext != ".nii.gz" && ext != ".nrrd") outFile += ".nii.gz"; mitk::IOUtil::Save(img, outFile); } delete handler; return EXIT_SUCCESS; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox index 1f9d80a938..d495c8bb24 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/documentation/UserManual/QmitkStreamlineTrackingViewUserManual.dox @@ -1,58 +1,59 @@ /** \page org_mitk_views_streamlinetracking Streamline Tractography This view enables streamline tractography on various input data. The corresponding command line application is named "MitkStreamlineTractography". Available sections: - \ref StrTrackUserManualInputData - \ref StrTrackUserManualInteractive - \ref StrTrackUserManualParameters - \ref StrTrackUserManualAdvancedParameters - \ref StrTrackUserManualReferences \section StrTrackUserManualInputData Input Data Input: \li One or multiple DTI images selected in the datamanager. \li One ODF image (e.g. obtained using MITK Q-ball reconstruction or MRtrix CSD). \li One peak image in MRtrix format (4D float image). Optional Input: \li Binary mask used to define the seed voxels. If no seed mask is specified, the whole image volume is seeded. \li Binary mask used to constrain the generated streamlines. Streamlines can not leave the mask area. \li Binary mask used to define stopping regions. Streamlines that enter the mask area are stopped immediately. \li Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. \li FA/GFA image used to determine streamline termination. If no image is specified, the FA/GFA image is automatically calculated from the input image. If multiple tensor images are used as input, it is recommended to provide such an image since the FA maps calculated from the individual input tensor images can not provide a suitable termination criterion. \section StrTrackUserManualInteractive Interactive Tractography Interactive tractography enables the dynamic placement of spherical seed regions simply by clicking into the image. Parameters are the number of seed points and the radius of the spherical seed region. The seed points are randomly distributed inside the sphere around the selected position. By clicking and holding the left mouse button while moving the mouse around the image, the fiber connections originating from the selected region can be explored dynamically. \section StrTrackUserManualParameters Parameters \li Mode: Toggle between deterministic and probabilistic tractography. Peak tracking only supports deterministic mode. The probabilistic method simply samples the output direction from the discrete probability ditribution provided by the discretized ODF. \li Seeds per voxel: If set to 1, the seed is defined as the voxel center. If > 1 the seeds are distributet randomly inside the voxel. \li Max. num. fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. -\li Cutoff: If the streamline reaches a position with an FA value or peak magnitude lower than the speciefied threshold, tracking is terminated. +\li Cutoff: If the streamline reaches a position with an FA value or peak magnitude lower than the speciefied threshold, tracking is terminated. Typical values are 0.2 for FA/GFA and 0.1 for CSD peaks. +\li ODF Cutoff: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. For MRtrix CSD fODF images, a typical value is 0.1. +\li Sharpen ODFs: If you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). This is not necessary (and not recommended) for CSD fODFs, since they are naturally much sharper. \section StrTrackUserManualAdvancedParameters Advanced Parameters \li Step Size: The algorithm proceeds along the streamline with a fixed stepsize. Default is 0.5*minSpacing. \li Angular threshold: Maximum angle between two successive steps (in degree). Default is 90° * step_size. For probabilistic tractography, candidate directions exceeding this threshold have probability 0, i.e. the respective ODF value is set to zero. The probabilities of the valid directions are normalized to sum to 1. \li Min. Tract Length: Shorter fibers are discarded. -\li Num. Samples per Step: Only relevant for probabilistic tractography. At each streamline position, the ODF is sampled n times and the direction with the highest probability is used. Increasing this number approximates the result to the deterministic tractography. \li f and g values to balance between FACT [1] and TEND [2,3] tracking (only for tensor based tractography). For further information please refer to [2,3] \li Flip directions: Internally flips progression directions. This might be necessary depending on the input data. \li Neighborhood Samples: Number of neighborhood samples that are used to determine the next fiber progression direction [4]. \li Compress Fibers: Whole brain tractograms obtained with a small step size can contain billions of points. The tractograms can be compressed by removing points that do not really contribute to the fiber shape, such as many points on a straight line. An error threshold (in mm) can be defined to specify which points should be removed and which not. \li Enable Trilinear Interpolation: By default the image values are interpolated. Keep in mind that in the noninterpolated case, the TEND term is only applied once per voxel. In the interpolated case the TEND term is applied at each integration step which results in much higher curvatures and has to be compensated by an according choice of f and g. \li Output Probability Map: No streamline are generated. Instead, the tractography outputs a probability map that indicates the probability of a fiber to reach a voxel from the selected seed region. For this measure to be sensible, the number of seeds per voxel needs to be rather large. This option does not work during interactive tractography. \li Enable Gray Matter Seeding: Seeds are onyl placed inside of the gray matter. Needs tissue label image. \section StrTrackUserManualReferences References [1] Mori et al. Annals Neurology 1999\n [2] Weinstein et al. Proceedings of IEEE Visualization 1999\n [3] Lazar et al. Human Brain Mapping 2003\n */ 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 6e470e2f33..2a2df6ef69 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,588 +1,592 @@ /*=================================================================== 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; QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : m_Controls(nullptr) , m_TrackingHandler(nullptr) { } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); m_Controls->m_FaImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_SeedImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_MaskImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_StopImageBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_TissueImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageBox->SetZeroEntryText("--"); m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_SeedImageBox->SetZeroEntryText("--"); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_StopImageBox->SetZeroEntryText("--"); m_Controls->m_TissueImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_TissueImageBox->SetZeroEntryText("--"); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_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()) ); m_FirstTensorProbRun = true; } UpdateGui(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); mitk::IRenderWindowPart* renderWindow = this->GetRenderWindowPart(); if ( m_Controls->m_InteractiveBox->isChecked() ) { QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractiveNode = mitk::DataNode::New(); QString name("InteractiveFib"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); 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; 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*/) { 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::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (m_TrackingHandler != nullptr) { delete m_TrackingHandler; m_TrackingHandler = nullptr; } } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer part, const QList& nodes ) { m_InputImageNodes.clear(); m_InputImages.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 { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr) { int dim = img->GetDimension(); unsigned int* dimensions = img->GetDimensions(); if (dim==4 && dimensions[3]%3==0) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); } } } } } UpdateGui(); } void QmitkStreamlineTrackingView::UpdateGui() { m_Controls->m_TensorImageLabel->setText("mandatory"); m_Controls->m_fBox->setVisible(false); m_Controls->m_fLabel->setVisible(false); m_Controls->m_gBox->setVisible(false); m_Controls->m_gLabel->setVisible(false); m_Controls->m_FaImageBox->setVisible(false); m_Controls->mFaImageLabel->setVisible(false); - m_Controls->m_numProbSamplesBox->setVisible(false); - m_Controls->m_numProbSamplesLabel->setVisible(false); + m_Controls->m_OdfCutoffBox->setVisible(false); + m_Controls->m_OdfCutoffLabel->setVisible(false); + m_Controls->m_SharpenOdfsBox->setVisible(false); if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_SeedGmBox->setVisible(true); else m_Controls->m_SeedGmBox->setVisible(false); if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.size()+" images selected"); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->m_InputData->setTitle("Input Data"); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_InteractiveBox->setEnabled(true); - if (m_Controls->m_ModeBox->currentIndex()==1) - { - m_Controls->m_numProbSamplesBox->setVisible(true); - m_Controls->m_numProbSamplesLabel->setVisible(true); - } - if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setVisible(true); m_Controls->m_fLabel->setVisible(true); m_Controls->m_gBox->setVisible(true); m_Controls->m_gLabel->setVisible(true); m_Controls->mFaImageLabel->setVisible(true); m_Controls->m_FaImageBox->setVisible(true); + +// if (m_Controls->m_ModeBox->currentIndex()==1) +// { +// m_Controls->m_OdfCutoffBox->setVisible(true); +// m_Controls->m_OdfCutoffLabel->setVisible(true); +// m_Controls->m_SharpenOdfsBox->setVisible(true); +// } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->mFaImageLabel->setVisible(true); m_Controls->m_FaImageBox->setVisible(true); - } - else - { + m_Controls->m_OdfCutoffBox->setVisible(true); + m_Controls->m_OdfCutoffLabel->setVisible(true); + m_Controls->m_SharpenOdfsBox->setVisible(true); } } else { m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_InteractiveBox->setEnabled(false); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_InputImages.empty()) return; if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; typedef mitk::ImageToItk CasterType; if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImages.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); return; } if (m_FirstTensorProbRun) { - QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. TEND parameters are ignored."); + 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)->SetNumProbSamples(m_Controls->m_numProbSamplesBox->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)->SetNumProbSamples(m_Controls->m_numProbSamplesBox->value()); + dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); + dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else { if (m_Controls->m_ModeBox->currentIndex()==1) { QMessageBox::information(nullptr, "Information", "Probabilstic tractography is not implemented for peak images."); 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(...) { 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); } typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if (m_Controls->m_InteractiveBox->isChecked()) { 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); } if (m_Controls->m_MaskImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_MaskImageBox->GetSelectedNode()->GetData()), mask); tracker->SetMaskImage(mask); } if (m_Controls->m_StopImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_StopImageBox->GetSelectedNode()->GetData()), mask); tracker->SetStoppingRegions(mask); } if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_TissueImageBox->GetSelectedNode()->GetData()), mask); tracker->SetTissueImage(mask); tracker->SetSeedOnlyGm(m_Controls->m_SeedGmBox->isChecked()); } tracker->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelBox->value()); tracker->SetStepSize(m_Controls->m_StepSizeBox->value()); //tracker->SetSamplingDistance(0.7); tracker->SetUseStopVotes(false); tracker->SetOnlyForwardSamples(false); tracker->SetAposterioriCurvCheck(false); tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); tracker->SetTrackingHandler(m_TrackingHandler); tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked() && !m_Controls->m_InteractiveBox->isChecked()); tracker->Update(); if (m_Controls->m_InteractiveBox->isChecked()) { vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("InteractiveFib"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } 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()); m_InteractiveNode->SetData(fib); } else if (!tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); if (fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); - warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter. Especially probabilistic tractography with a low number of samples per step tends to produce very curvy tracts that are then rejected.\n- Number of samples in probabilistic tractography is too low. See angular threshold bullet.\n- A small step sizes also means many steps to go wrong. See angular threshold bullet."); + 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(); 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 = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_InputImageNodes.at(0)->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_InputImageNodes.at(0)); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingViewControls.ui index e6bbac6688..61f288feff 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,787 +1,800 @@ QmitkStreamlineTrackingViewControls 0 0 406 950 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 - 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: Input Image. ODF, tensor and peak images are currently supported. <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Seed points are only placed inside the mask image. If no seed mask is selected, the whole image is seeded. QComboBox::AdjustToMinimumContentsLength - Stop Image: Input Image. ODF, tensor and peak images are currently supported. Tensor/Peak/ODF Image: FA/GFA Image: Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. QComboBox::AdjustToMinimumContentsLength - Tissue Image: Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. QComboBox::AdjustToMinimumContentsLength - false Start Tractography 0 0 Parameters Mode: Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. Deterministic Probabilistic Seeds per Voxel: Number of seed points placed in each voxel. 1 9999999 Max. num. fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. 999999999 Cutoff: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 - + Qt::Horizontal QSizePolicy::Fixed 200 0 + + + + + + + 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 + + + Advanced Paramaters Step Size: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Angular Threshold: Default: 90° * step_size -1 90 1 -1 Min. Tract Length: Shorter fibers are discarded. Minimum fiber length (in mm) 1 999.000000000000000 1.000000000000000 20.000000000000000 - + f parameter of tensor tractography. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). f: - + f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 1.000000000000000 - + g: - + f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 0.000000000000000 - + Flip directions: - + QFrame::NoFrame QFrame::Raised 0 0 0 0 Internally flips progression directions. This might be necessary depending on the input data. x Internally flips progression directions. This might be necessary depending on the input data. y Internally flips progression directions. This might be necessary depending on the input data. z - + Neighborhood Samples: - + Number of neighborhood samples that are used to determine the next fiber progression direction. 50 - + Resample fibers using the specified error constraint. Compress Fibers true - + Qt::StrongFocus Lossy fiber compression. Recommended for large tractograms. Maximum error in mm. 3 10.000000000000000 0.010000000000000 0.100000000000000 - + If false, nearest neighbor interpolation is used. Enable Trilinear Interpolation true - + Requires tissue image. Enable Gray Matter Seeding false - + Output probability map instead of tractogram. Output Probability Map false - - - - 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!). - - - Num. Samples per Step: - - - - - - - At each integration step the ODF is sampled n times and the most probable direction is used. The higher the number, the closer the behaviour will be to the deterministic ODF tractography. - - - 1 - - - 9999999 - - - 10 - - - Qt::Vertical QSizePolicy::Expanding 20 220 0 0 Interactive Tractography Radius: false Dynamically pick seed location by click into image. Enable Interactive Tractography Num. Seeds: Number of seed points normally distributed around selected position. 1 9999999 50 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 QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h