diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx index 1d43ded4d5..8b65d51254 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx @@ -1,864 +1,864 @@ /*=================================================================== 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 _itk_TensorReconstructionWithEigenvalueCorrectionFilter_txx_ #define _itk_TensorReconstructionWithEigenvalueCorrectionFilter_txx_ #endif #include "itkImageRegionConstIterator.h" #include #include "itkImageFileWriter.h" #include "itkImage.h" #include "itkImageRegionIterator.h" #include namespace itk { template TensorReconstructionWithEigenvalueCorrectionFilter ::TensorReconstructionWithEigenvalueCorrectionFilter() { m_B0Threshold = 50.0; } template void TensorReconstructionWithEigenvalueCorrectionFilter ::GenerateData () { typename GradientImagesType::Pointer input_image = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); typename itk::ImageDuplicator::Pointer duplicator = itk::ImageDuplicator::New(); duplicator->SetInputImage(input_image); duplicator->Update(); m_GradientImagePointer = duplicator->GetOutput(); typename GradientImagesType::SizeType size = m_GradientImagePointer->GetLargestPossibleRegion().GetSize(); // number of volumes int nof = m_GradientDirectionContainer->Size(); // determine the number of b-zero values int numberb0=0; for(int i=0; i vec = m_GradientDirectionContainer->ElementAt(i); float bval = vec.magnitude(); bval = bval*bval*m_BValue; if(bval<100) numberb0++; } // Matrix to store all diffusion encoding gradients vnl_matrix directions(nof-numberb0,3); m_B0Mask.set_size(nof); int cnt=0; for(int i=0; i vec = m_GradientDirectionContainer->ElementAt(i); float bval = vec.magnitude(); bval = bval*bval*m_BValue; if(bval<100) { // the diffusion encoding gradient is approximately zero, wo we are dealing with a non-diffusion weighted volume m_B0Mask[i]=1; } else { // dealing with a diffusion weighted volume m_B0Mask[i]=0; // set the diffusion encoding gradient to the directions matrix directions[cnt][0] = vec[0]; directions[cnt][1] = vec[1]; directions[cnt][2] = vec[2]; cnt++; } } // looking for maximal norm among gradients. // The norm is calculated with use of spectral radius theorem- based on determination of eigenvalue. vnl_matrix dirsTimesDirsTrans = directions*directions.transpose(); vnl_vector< double> diagonal(nof-numberb0); vnl_vector< double> b_vec(nof-numberb0); vnl_vector< double> temporary(3); for (int i=0;i H(nof-numberb0, 6); vnl_matrix H_org(nof-numberb0, 6); vnl_vector pre_tensor(9); //H is matrix that contains covariances for directions. It is stored twice because its original value is needed later // while H is changed int etbt[6] = { 0, 4, 8, 1, 5, 2 };// tensor order for (int i = 0; i < nof-numberb0; i++) { for (int j = 0; j < 3; j++) { temporary[j] = -directions[i][j]; } for (int j = 0; j < 3; j++) { for (int k = 0; k < 3; k++) { pre_tensor[k + 3 * j] = temporary[k] * directions[i][j]; } } for (int j = 0; j < 6; j++) { H[i][j] = pre_tensor[etbt[j]]; } for (int j = 0; j < 3; j++) { H[i][3 + j] *= 2.0; } } H_org=H; // calculation of inverse matrix by means of pseudoinverse vnl_matrix inputtopseudoinverse=H.transpose()*H; vnl_symmetric_eigensystem eig( inputtopseudoinverse); vnl_matrix pseudoInverse = eig.pinverse()*H.transpose(); typedef itk::Image MaskImageType; MaskImageType::Pointer mask = MaskImageType::New(); mask->SetRegions(m_GradientImagePointer->GetLargestPossibleRegion().GetSize()); mask->SetSpacing(m_GradientImagePointer->GetSpacing()); mask->SetOrigin(m_GradientImagePointer->GetOrigin()); mask->SetDirection( m_GradientImagePointer->GetDirection() ); // Set the image direction mask->SetLargestPossibleRegion( m_GradientImagePointer->GetLargestPossibleRegion() ); mask->SetBufferedRegion( m_GradientImagePointer->GetLargestPossibleRegion() ); mask->SetRequestedRegion( m_GradientImagePointer->GetLargestPossibleRegion() ); mask->Allocate(); // Image thresholding: For every voxel mean B0 image is calculated and then voxels of mean B0 less than the // treshold on the B0 image proviced by the userare excluded from the dataset with use of defined mask image. // 1 in mask voxel means that B0 > assumed treshold. int mask_cnt=0; #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif - for(itk::SizeValueType x=0;x ix = {{(itk::IndexValueType)x,(itk::IndexValueType)y,(itk::IndexValueType)z}}; GradientVectorType pixel = m_GradientImagePointer->GetPixel(ix); for (int i=0;i m_B0Threshold) { #pragma omp critical { mask->SetPixel(ix, 1); mask_cnt++; } } else { #pragma omp critical mask->SetPixel(ix, 0); } } #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif - for (itk::SizeValueType x=0;x org_vec(nof); itk::Index<3> ix = {{(itk::IndexValueType)x,(itk::IndexValueType)y,(itk::IndexValueType)z}}; double mask_val = mask->GetPixel(ix); GradientVectorType pixel2 = m_GradientImagePointer->GetPixel(ix); for (int i=0;i0) { for( int f=0;fSetPixel(ix, pixel2); } } typename TensorImageType::Pointer tensorImg; tensorImg = TensorImageType::New(); tensorImg->SetSpacing( m_GradientImagePointer->GetSpacing() ); // Set the image spacing tensorImg->SetOrigin( m_GradientImagePointer->GetOrigin() ); // Set the image origin tensorImg->SetDirection( m_GradientImagePointer->GetDirection() ); // Set the image direction tensorImg->SetLargestPossibleRegion( m_GradientImagePointer->GetLargestPossibleRegion() ); tensorImg->SetBufferedRegion( m_GradientImagePointer->GetLargestPossibleRegion() ); tensorImg->SetRequestedRegion( m_GradientImagePointer->GetLargestPossibleRegion() ); tensorImg->Allocate(); //Declaration of vectors that contains too high or too low atenuation for each gradient. Attenuation is only calculated for //non B0 images so nof-numberb0. vnl_vector< double> pixel_max(nof-numberb0); vnl_vector< double> pixel_min(nof-numberb0); // to high and to low attenuation is calculated with use of highest allowed =5 and lowest allowed =0.01 diffusion coefficient for (int i=0;i outputIterator(tensorImg, tensorImg->GetLargestPossibleRegion()); outputIterator.GoToBegin(); while(!outputIterator.IsAtEnd()) { TensorPixelType tens = outputIterator.Get(); tens/= 1000.0; outputIterator.Set(tens); ++outputIterator; } this->SetNthOutput(0, tensorImg); } template void TensorReconstructionWithEigenvalueCorrectionFilter ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) { itkExceptionMacro( << "Cannot call both methods:" << "AddGradientImage and SetGradientImage. Please call only one of them."); } this->m_GradientDirectionContainer = gradientDirection; unsigned int numImages = gradientDirection->Size(); this->m_NumberOfBaselineImages = 0; this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; // ensure that the gradient image we received has as many components as // the number of gradient directions if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + this->m_NumberOfGradientDirections ) { itkExceptionMacro( << this->m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages << "baselines = " << this->m_NumberOfGradientDirections + this->m_NumberOfBaselineImages << " directions specified but image has " << gradientImage->GetVectorLength() << " components."); } this->ProcessObject::SetNthInput( 0, const_cast< GradientImagesType* >(gradientImage) ); m_GradientImageTypeEnumeration = GradientIsInASingleImage; } template double TensorReconstructionWithEigenvalueCorrectionFilter ::CheckNeighbours(int x, int y, int z,int f, itk::Size<3> size, itk::Image::Pointer mask, typename GradientImagesType::Pointer corrected_diffusion_temp) { // method is used for finding a new value for the voxel with use of its 27 neighborhood. To perform such a smoothing correct voxels are // counted an arithmetical mean is calculated and stored as a new value for the voxel. If there is no proper neigborhood voxel is turned // to the value of 0. // Definition of neighbourhood avoiding crossing the image boundaries int x_max=size[0]-1; int y_max=size[1]-1; int z_max=size[2]-1; double back_x=std::max(0,x-1); double back_y=std::max(0,y-1); double back_z=std::max(0,z-1); double forth_x=std::min((x+1),x_max); double forth_y=std::min((y+1),y_max); double forth_z=std::min((z+1),z_max); double tempsum=0; double temp_number=0; for(int i=back_x; i<=forth_x; i++) { for (int j=back_y; j<=forth_y; j++) { for (int k=back_z; k<=forth_z; k++) { itk::Index<3> ix = {{i,j,k}}; GradientVectorType p = corrected_diffusion_temp->GetPixel(ix); if (p[f] > 0.0 )// taking only positive values and counting them { if(!(i==x && j==y && k== z)) { tempsum=tempsum+p[f]; temp_number++; } } } } } //getting back to the original position of the voxel itk::Index<3> ix = {{x,y,z}}; if (temp_number <= 0.0) { tempsum=0; #pragma omp critical mask->SetPixel(ix,0); } else { tempsum=tempsum/temp_number; } return tempsum;// smoothed value of voxel } template void TensorReconstructionWithEigenvalueCorrectionFilter ::CalculateAttenuation(vnl_vector org_data,vnl_vector &atten,int nof, int numberb0) { double mean_b=0.0; for (int i=0;i0) mean_b=mean_b+org_data[i]; mean_b=mean_b/numberb0; int cnt=0; for (int i=0;i double TensorReconstructionWithEigenvalueCorrectionFilter ::CheckNegatives ( itk::Size<3> size, itk::Image::Pointer mask, typename itk::Image< itk::DiffusionTensor3D, 3 >::Pointer tensorImg ) { // The method was created to simplif the flow of negative eigenvalue correction process. The method itself just return the number // of voxels (tensors) with negative eigenvalues. Then if the voxel was previously bad ( mask=2 ) but it is not bad anymore mask is //changed to 1. // declaration of important structures and variables double badvoxels=0; #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif - for (itk::SizeValueType x=0;x ix = {{(itk::IndexValueType)x,(itk::IndexValueType)y,(itk::IndexValueType)z}}; pixel = mask->GetPixel(ix); // but only if previously marked as bad one-negative eigen value if(pixel > 1) { #pragma omp critical { itk::DiffusionTensor3D::EigenValuesArrayType eigenvalues; itk::DiffusionTensor3D::EigenVectorsMatrixType eigenvectors; itk::DiffusionTensor3D ten = tensorImg->GetPixel(ix); ten.ComputeEigenAnalysis(eigenvalues, eigenvectors); //comparison to 0.01 instead of 0 was proposed by O.Pasternak if( eigenvalues[0]>0.01 && eigenvalues[1]>0.01 && eigenvalues[2]>0.01) mask->SetPixel(ix,1); else badvoxels++; } } } return badvoxels; } template void TensorReconstructionWithEigenvalueCorrectionFilter ::CorrectDiffusionImage(int nof, int numberb0, itk::Size<3> size, typename GradientImagesType::Pointer corrected_diffusion,itk::Image::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min) { // in this method the voxels that has tensor negative eigenvalues are smoothed. Smoothing is done on DWI image.For the voxel //detected as bad one, B0 image is smoothed obligatory. All other gradient images are smoothed only when value of attenuation //is out of declared bounds for too high or too low attenuation. // declaration of important variables #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif - for (itk::SizeValueType z=0;z org_data(nof); vnl_vector atten(nof-numberb0); double cnt_atten=0; itk::Index<3> ix = {{(itk::IndexValueType)x, (itk::IndexValueType)y, (itk::IndexValueType)z}}; if(mask->GetPixel(ix) > 1.0) { GradientVectorType pt = corrected_diffusion->GetPixel(ix); for (int i=0;i0) { mean_b=mean_b+org_data[i]; } mean_b=mean_b/numberb0; int cnt=0; for (int i=0;i pixel_max[cnt_atten]) { int x_max=size[0]-1; int y_max=size[1]-1; int z_max=size[2]-1; double back_x=std::max(0,(int)x-1); double back_y=std::max(0,(int)y-1); double back_z=std::max(0,(int)z-1); double forth_x=std::min(((int)x+1),x_max); double forth_y=std::min(((int)y+1),y_max); double forth_z=std::min(((int)z+1),z_max); double tempsum=0; double temp_number=0; - for(unsigned int i=back_x; i<=forth_x; i++) - for (unsigned int j=back_y; j<=forth_y; j++) - for (unsigned int k=back_z; k<=forth_z; k++) + for(int i=back_x; i<=forth_x; i++) + for (int j=back_y; j<=forth_y; j++) + for (int k=back_z; k<=forth_z; k++) { itk::Index<3> ix = {{i,j,k}}; GradientVectorType p = corrected_diffusion->GetPixel(ix); if(p[f] > 0.0 && !(i==x && j==y && k== z)) { tempsum=tempsum+p[f]; temp_number++; } } //getting back to the original position of the voxel itk::Index<3> ix = {{(itk::IndexValueType)x,(itk::IndexValueType)y,(itk::IndexValueType)z}}; if (temp_number <= 0.0) { tempsum=0; #pragma omp critical mask->SetPixel(ix,0); } else tempsum=tempsum/temp_number; org_data[f] = tempsum; } cnt_atten++; } //smoothing B0 if(m_B0Mask[f]==1) { int x_max=size[0] - 1; int y_max=size[1] - 1; int z_max=size[2] - 1; double back_x=std::max(0,(int)x-1); double back_y=std::max(0,(int)y-1); double back_z=std::max(0,(int)z-1); double forth_x=std::min(((int)x+1),x_max); double forth_y=std::min(((int)y+1),y_max); double forth_z=std::min(((int)z+1),z_max); double tempsum=0; double temp_number=0; - for(unsigned int i=back_x; i<=forth_x; i++) - for (unsigned int j=back_y; j<=forth_y; j++) - for (unsigned int k=back_z; k<=forth_z; k++) + for(int i=back_x; i<=forth_x; i++) + for (int j=back_y; j<=forth_y; j++) + for (int k=back_z; k<=forth_z; k++) { itk::Index<3> ix = {{i,j,k}}; GradientVectorType p = corrected_diffusion->GetPixel(ix); //double test= p[f]; if (p[f] > 0.0 )// taking only positive values and counting them { if(!(i==x && j==y && k== z)) { tempsum=tempsum+p[f]; temp_number++; } } } //getting back to the original position of the voxel itk::Index<3> ix = {{(itk::IndexValueType)x,(itk::IndexValueType)y,(itk::IndexValueType)z}}; if (temp_number <= 0.0) { tempsum=0; #pragma omp critical mask->SetPixel(ix,0); } else { tempsum=tempsum/temp_number; } org_data[f] = tempsum; } } for (int i=0;iSetPixel(ix, pt); } else { GradientVectorType pt = corrected_diffusion->GetPixel(ix); #pragma omp critical corrected_diffusion->SetPixel(ix, pt); } } } template void TensorReconstructionWithEigenvalueCorrectionFilter ::GenerateTensorImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,double , typename itk::Image< itk::DiffusionTensor3D, 3 >::Pointer tensorImg) { // in this method the whole tensor image is updated with a tensors for defined voxels ( defined by a value of mask); #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif - for (itk::SizeValueType x=0;x ix; vnl_vector org_data(nof); vnl_vector atten(nof-numberb0); vnl_vector tensor(6); itk::DiffusionTensor3D ten; double mask_val=0; ix[0] = x; ix[1] = y; ix[2] = z; mask_val= mask->GetPixel(ix); //Tensors are calculated only for voxels above theshold for B0 image. if( mask_val > 0.0 ) { // calculation of attenuation with use of gradient image and and mean B0 image GradientVectorType pt = corrected_diffusion->GetPixel(ix); for (int i=0;i0) { mean_b=mean_b+org_data[i]; } } mean_b=mean_b/numberb0; int cnt=0; for (int i=0;iSetPixel(ix, ten); } // for voxels with mask value 0 - tensor is simply 0 ( outside brain value) else if (mask_val < 1.0) { ten(0,0) = 0; ten(0,1) = 0; ten(0,2) = 0; ten(1,1) = 0; ten(1,2) = 0; ten(2,2) = 0; #pragma omp critical tensorImg->SetPixel(ix, ten); } } }// end of Generate Tensor template void TensorReconstructionWithEigenvalueCorrectionFilter ::TurnMask( itk::Size<3> size, itk::Image::Pointer mask, double previous_mask, double set_mask) { // The method changes voxels in the mask that poses a certain value with other value. itk::Index<3> ix; double temp_mask_value=0; #ifdef WIN32 #pragma omp parallel for #else #pragma omp parallel for collapse(3) #endif - for(itk::SizeValueType x=0;xGetPixel(ix); if(temp_mask_value>previous_mask) { #pragma omp critical mask->SetPixel(ix,set_mask); } } } } // end of namespace 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 0b7eb085af..e19a991a1f 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,61 +1,74 @@ /** \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: +Select the data you want to track on in the datamanager. Supported file types are: \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). +\li One ODF image, e.g. obtained using MITK Q-ball reconstruction or MRtrix CSD (tractography similar to [6]). +\li One peak image (4D float image). +\li One raw diffusion-weighted image for machine learning based tractography [1]. 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. +\li Tractography Forest: Needed for machine learning based tractography [1]. \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. - +Interactive tractography enables the dynamic placement of spherical seed regions simply by clicking into the image (similar to [5]). 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. When "Update on Parameter Change" is checked, each parameter change causes an instant retracking with the new parameters. This enables an intuitive exploration of the effcts that the individual parameters have on the resulting tractogram. \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. 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 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 f and g values to balance between FACT [2] and TEND [3,4] 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 StrTrackUserManualAdvancedParameters Neighbourhood Sampling (for details see [1]) +\li Neighborhood Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. +\li Sampling Distance: Distance of the sampling positions from the current streamline position (in voxels). +\li Use Only Frontal Samples: Only neighborhood samples in front of the current streamline position are considered. +\li Use Stop-Votes: If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. + +\section StrTrackUserManualAdvancedParameters Output and Postprocessing +\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 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. + \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 +[1] Neher, Peter F., Marc-Alexandre Côté, Jean-Christophe Houde, Maxime Descoteaux, and Klaus H. Maier-Hein. “Fiber Tractography Using Machine Learning.” NeuroImage. Accessed July 19, 2017. doi:10.1016/j.neuroimage.2017.07.028.\n +[2] Mori, Susumu, Walter E. Kaufmann, Godfrey D. Pearlson, Barbara J. Crain, Bram Stieltjes, Meiyappan Solaiyappan, and Peter C. M. Van Zijl. “In Vivo Visualization of Human Neural Pathways by Magnetic Resonance Imaging.” Annals of Neurology 47 (2000): 412–414.\n +[3] Weinstein, David, Gordon Kindlmann, and Eric Lundberg. “Tensorlines: Advection-Diffusion Based Propagation through Diffusion Tensor Fields.” In Proceedings of the Conference on Visualization’99: Celebrating Ten Years, 249–253, n.d.\n +[4] Lazar, Mariana, David M. Weinstein, Jay S. Tsuruda, Khader M. Hasan, Konstantinos Arfanakis, M. Elizabeth Meyerand, Benham Badie, et al. “White Matter Tractography Using Diffusion Tensor Deflection.” Human Brain Mapping 18, no. 4 (2003): 306–321.\n +[5] Chamberland, M., K. Whittingstall, D. Fortin, D. Mathieu, and M. Descoteaux. “Real-Time Multi-Peak Tractography for Instantaneous Connectivity Display.” Front Neuroinform 8 (2014): 59. doi:10.3389/fninf.2014.00059.\n +[6] Tournier, J-Donald, Fernando Calamante, and Alan Connelly. “MRtrix: Diffusion Tractography in Crossing Fiber Regions.” International Journal of Imaging Systems and Technology 22, no. 1 (March 2012): 53–66. doi:10.1002/ima.22005. + */ 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 e70f90259b..d83356dee7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tractography/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,747 +1,747 @@ /*=================================================================== 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()); m_Controls->m_ForestBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isImagePredicate = mitk::TNodePredicateDataType::New(); mitk::TNodePredicateDataType::Pointer isTractographyForest = mitk::TNodePredicateDataType::New(); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); mitk::NodePredicateNot::Pointer isNotBinaryPredicate = mitk::NodePredicateNot::New( isBinaryPredicate ); mitk::NodePredicateAnd::Pointer isNotABinaryImagePredicate = mitk::NodePredicateAnd::New( isImagePredicate, isNotBinaryPredicate ); mitk::NodePredicateDimension::Pointer dimensionPredicate = mitk::NodePredicateDimension::New(3); m_Controls->m_ForestBox->SetPredicate(isTractographyForest); m_Controls->m_FaImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_FaImageBox->SetZeroEntryText("--"); m_Controls->m_SeedImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_SeedImageBox->SetZeroEntryText("--"); m_Controls->m_MaskImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_MaskImageBox->SetZeroEntryText("--"); m_Controls->m_StopImageBox->SetPredicate( mitk::NodePredicateAnd::New(isBinaryPredicate, dimensionPredicate) ); m_Controls->m_StopImageBox->SetZeroEntryText("--"); m_Controls->m_TissueImageBox->SetPredicate( mitk::NodePredicateAnd::New(isNotABinaryImagePredicate, dimensionPredicate) ); m_Controls->m_TissueImageBox->SetZeroEntryText("--"); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_InteractiveBox, SIGNAL(stateChanged(int)), this, SLOT(ToggleInteractive()) ); connect( m_Controls->m_TissueImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(UpdateGui()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(DeleteTrackingHandler()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OutputStyleSwitched()) ); connect( m_Controls->m_ModeBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MaskImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FaImageBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(ForestSwitched()) ); connect( m_Controls->m_ForestBox, SIGNAL(currentIndexChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedsPerVoxelBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumFibersBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_ScalarThresholdBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_OdfCutoffBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StepSizeBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SamplingDistanceBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_AngularThresholdBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_MinTractLengthBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_fBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_gBox, SIGNAL(valueChanged(double)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_NumSamplesBox, SIGNAL(valueChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedRadiusBox, SIGNAL(valueChanged(double)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_NumSeedsBox, SIGNAL(valueChanged(int)), this, SLOT(InteractiveSeedChanged()) ); connect( m_Controls->m_OutputProbMap, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SharpenOdfsBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_InterpolationBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_SeedGmBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipXBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipYBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FlipZBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_FrontalSamplesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); connect( m_Controls->m_StopVotesBox, SIGNAL(stateChanged(int)), this, SLOT(OnParameterChanged()) ); m_FirstTensorProbRun = true; } UpdateGui(); } void QmitkStreamlineTrackingView::InteractiveSeedChanged(bool posChanged) { if (!posChanged && (!m_Controls->m_InteractiveBox->isChecked() || !m_Controls->m_ParamUpdateBox->isChecked())) return; std::srand(std::time(0)); m_SeedPoints.clear(); itk::Point world_pos = this->GetRenderWindowPart()->GetSelectedPosition(); m_SeedPoints.push_back(world_pos); float radius = m_Controls->m_SeedRadiusBox->value(); int num = m_Controls->m_NumSeedsBox->value(); mitk::PointSet::Pointer pointset = mitk::PointSet::New(); pointset->InsertPoint(0, world_pos); m_InteractivePointSetNode->SetProperty("pointsize", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetProperty("point 2D size", mitk::FloatProperty::New(radius*2)); m_InteractivePointSetNode->SetData(pointset); for (int i=1; i p; p[0] = rand()%1000-500; p[1] = rand()%1000-500; p[2] = rand()%1000-500; p.Normalize(); p *= radius; m_SeedPoints.push_back(world_pos+p); } DoFiberTracking(); } void QmitkStreamlineTrackingView::OnParameterChanged() { if (m_Controls->m_InteractiveBox->isChecked() && m_Controls->m_ParamUpdateBox->isChecked()) DoFiberTracking(); } void QmitkStreamlineTrackingView::ToggleInteractive() { UpdateGui(); mitk::IRenderWindowPart* renderWindow = this->GetRenderWindowPart(); m_Controls->m_SeedsPerVoxelBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedsPerVoxelLabel->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedGmBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_SeedImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_6->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_TissueImageBox->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->label_10->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); if ( m_Controls->m_InteractiveBox->isChecked() ) { QApplication::setOverrideCursor(Qt::PointingHandCursor); QApplication::processEvents(); m_InteractivePointSetNode = mitk::DataNode::New(); m_InteractivePointSetNode->SetProperty("color", mitk::ColorProperty::New(1,1,1)); m_InteractivePointSetNode->SetName("InteractiveSeedRegion"); mitk::PointSetShapeProperty::Pointer shape_prop = mitk::PointSetShapeProperty::New(); shape_prop->SetValue(mitk::PointSetShapeProperty::PointSetShape::CIRCLE); m_InteractivePointSetNode->SetProperty("Pointset.2D.shape", shape_prop); GetDataStorage()->Add(m_InteractivePointSetNode); if (renderWindow) { { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("axial"))->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); m_SliceObserverTag1 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); } { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("sagittal"))->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); m_SliceObserverTag2 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); } { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("coronal"))->GetSliceNavigationController(); itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); command->SetCallbackFunction( this, &QmitkStreamlineTrackingView::OnSliceChanged ); m_SliceObserverTag3 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(nullptr, 0), command ); } } } else { QApplication::restoreOverrideCursor(); QApplication::processEvents(); m_InteractiveNode = nullptr; m_InteractivePointSetNode = nullptr; if (renderWindow) { mitk::SliceNavigationController* slicer = renderWindow->GetQmitkRenderWindow(QString("axial"))->GetSliceNavigationController(); slicer->RemoveObserver(m_SliceObserverTag1); slicer = renderWindow->GetQmitkRenderWindow(QString("sagittal"))->GetSliceNavigationController(); slicer->RemoveObserver(m_SliceObserverTag2); slicer = renderWindow->GetQmitkRenderWindow(QString("coronal"))->GetSliceNavigationController(); slicer->RemoveObserver(m_SliceObserverTag3); } } } void QmitkStreamlineTrackingView::OnSliceChanged(const itk::EventObject& /*e*/) { InteractiveSeedChanged(true); } void QmitkStreamlineTrackingView::SetFocus() { } void QmitkStreamlineTrackingView::DeleteTrackingHandler() { if (m_TrackingHandler != nullptr) { delete m_TrackingHandler; m_TrackingHandler = nullptr; } } void QmitkStreamlineTrackingView::ForestSwitched() { DeleteTrackingHandler(); } void QmitkStreamlineTrackingView::OutputStyleSwitched() { if (m_InteractiveNode.IsNotNull()) GetDataStorage()->Remove(m_InteractiveNode); m_InteractiveNode = nullptr; } void QmitkStreamlineTrackingView::OnSelectionChanged( berry::IWorkbenchPart::Pointer , const QList& nodes ) { m_InputImageNodes.clear(); m_InputImages.clear(); m_AdditionalInputImages.clear(); DeleteTrackingHandler(); for( auto node : nodes ) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); } else if ( dynamic_cast(node->GetData()) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData())) ) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); } else { mitk::Image* img = dynamic_cast(node->GetData()); if (img!=nullptr) { int dim = img->GetDimension(); unsigned int* dimensions = img->GetDimensions(); if (dim==4 && dimensions[3]%3==0) { m_InputImageNodes.push_back(node); m_InputImages.push_back(dynamic_cast(node->GetData())); } else if (dim==3) { m_AdditionalInputImages.push_back(dynamic_cast(node->GetData())); } } } } } UpdateGui(); OnParameterChanged(); } void QmitkStreamlineTrackingView::UpdateGui() { m_Controls->m_TensorImageLabel->setText("mandatory"); m_Controls->m_fBox->setEnabled(false); m_Controls->m_fLabel->setEnabled(false); m_Controls->m_gBox->setEnabled(false); m_Controls->m_gLabel->setEnabled(false); m_Controls->m_FaImageBox->setEnabled(false); m_Controls->mFaImageLabel->setEnabled(false); m_Controls->m_OdfCutoffBox->setEnabled(false); m_Controls->m_OdfCutoffLabel->setEnabled(false); m_Controls->m_SharpenOdfsBox->setEnabled(false); m_Controls->m_ForestBox->setEnabled(false); m_Controls->m_ForestLabel->setEnabled(false); if (m_Controls->m_TissueImageBox->GetSelectedNode().IsNotNull()) m_Controls->m_SeedGmBox->setEnabled(true); else m_Controls->m_SeedGmBox->setEnabled(false); if(!m_InputImageNodes.empty()) { if (m_InputImageNodes.size()>1) m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.size()+" images selected"); else m_Controls->m_TensorImageLabel->setText(m_InputImageNodes.at(0)->GetName().c_str()); m_Controls->m_InputData->setTitle("Input Data"); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); m_Controls->m_ScalarThresholdBox->setEnabled(true); m_Controls->m_FaThresholdLabel->setEnabled(true); if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->m_fBox->setEnabled(true); m_Controls->m_fLabel->setEnabled(true); m_Controls->m_gBox->setEnabled(true); m_Controls->m_gLabel->setEnabled(true); m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { m_Controls->mFaImageLabel->setEnabled(true); m_Controls->m_FaImageBox->setEnabled(true); m_Controls->m_OdfCutoffBox->setEnabled(true); m_Controls->m_OdfCutoffLabel->setEnabled(true); m_Controls->m_SharpenOdfsBox->setEnabled(true); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { m_Controls->m_ForestBox->setEnabled(true); m_Controls->m_ForestLabel->setEnabled(true); m_Controls->m_ScalarThresholdBox->setEnabled(false); m_Controls->m_FaThresholdLabel->setEnabled(false); } } else { m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->commandLinkButton->setEnabled(!m_Controls->m_InteractiveBox->isChecked()); } } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_InputImages.empty()) return; typedef itk::StreamlineTrackingFilter TrackerType; TrackerType::Pointer tracker = TrackerType::New(); if( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; typedef mitk::ImageToItk CasterType; if (m_Controls->m_ModeBox->currentIndex()==1) { if (m_InputImages.size()>1) { QMessageBox::information(nullptr, "Information", "Probabilistic tensor tractography is only implemented for single-tensor mode!"); return; } if (m_FirstTensorProbRun) { QMessageBox::information(nullptr, "Information", "Internally calculating ODF from tensor image and performing probabilistic ODF tractography. ODFs are sharpened (min-max normalized and raised to the power of 4). TEND parameters are ignored."); m_FirstTensorProbRun = false; } if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); TensorImageType::Pointer itkImg = TensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); typedef itk::TensorImageToQBallImageFilter< float, float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkImg ); filter->Update(); dynamic_cast(m_TrackingHandler)->SetOdfImage(filter->GetOutput()); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(0); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(true); } else { if (m_TrackingHandler==nullptr) { m_TrackingHandler = new mitk::TrackingHandlerTensor(); for (int i=0; i<(int)m_InputImages.size(); i++) { TensorImageType::Pointer itkImg = TensorImageType::New(); mitk::CastToItkImage(m_InputImages.at(i), itkImg); dynamic_cast(m_TrackingHandler)->AddTensorImage(itkImg); } if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetFaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetFaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetF((float)m_Controls->m_fBox->value()); dynamic_cast(m_TrackingHandler)->SetG((float)m_Controls->m_gBox->value()); } } else if ( dynamic_cast(m_InputImageNodes.at(0)->GetData()) ) { if (m_TrackingHandler==nullptr) { typedef mitk::ImageToItk< mitk::TrackingHandlerOdf::ItkOdfImageType > CasterType; m_TrackingHandler = new mitk::TrackingHandlerOdf(); mitk::TrackingHandlerOdf::ItkOdfImageType::Pointer itkImg = mitk::TrackingHandlerOdf::ItkOdfImageType::New(); mitk::CastToItkImage(m_InputImages.at(0), itkImg); dynamic_cast(m_TrackingHandler)->SetOdfImage(itkImg); if (m_Controls->m_FaImageBox->GetSelectedNode().IsNotNull()) { ItkFloatImageType::Pointer itkImg = ItkFloatImageType::New(); mitk::CastToItkImage(dynamic_cast(m_Controls->m_FaImageBox->GetSelectedNode()->GetData()), itkImg); dynamic_cast(m_TrackingHandler)->SetGfaImage(itkImg); } } dynamic_cast(m_TrackingHandler)->SetGfaThreshold(m_Controls->m_ScalarThresholdBox->value()); dynamic_cast(m_TrackingHandler)->SetOdfThreshold(m_Controls->m_OdfCutoffBox->value()); dynamic_cast(m_TrackingHandler)->SetSharpenOdfs(m_Controls->m_SharpenOdfsBox->isChecked()); } else if ( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_InputImageNodes.at(0)->GetData())) ) { if ( m_Controls->m_ForestBox->GetSelectedNode().IsNull() ) { - QMessageBox::information(nullptr, "Information", "Not random forest for machine-learning based tractography selected."); + QMessageBox::information(nullptr, "Information", "Not random forest for machine learning based tractography selected."); return; } if (m_TrackingHandler==nullptr) { mitk::TractographyForest::Pointer forest = dynamic_cast(m_Controls->m_ForestBox->GetSelectedNode()->GetData()); mitk::Image::Pointer dwi = dynamic_cast(m_InputImageNodes.at(0)->GetData()); std::vector< std::vector< ItkFloatImageType::Pointer > > additionalFeatureImages; additionalFeatureImages.push_back(std::vector< ItkFloatImageType::Pointer >()); for (auto img : m_AdditionalInputImages) { ItkFloatImageType::Pointer itkimg = ItkFloatImageType::New(); mitk::CastToItkImage(img, itkimg); additionalFeatureImages.at(0).push_back(itkimg); } bool forest_valid = false; if (forest->GetNumFeatures()>=100) { int num_previous_directions = (forest->GetNumFeatures() - (100 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 100>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } else { int num_previous_directions = (forest->GetNumFeatures() - (28 + additionalFeatureImages.at(0).size()))/3; m_TrackingHandler = new mitk::TrackingHandlerRandomForest<6, 28>(); dynamic_cast*>(m_TrackingHandler)->AddDwi(dwi); dynamic_cast*>(m_TrackingHandler)->SetAdditionalFeatureImages(additionalFeatureImages); dynamic_cast*>(m_TrackingHandler)->SetForest(forest); dynamic_cast*>(m_TrackingHandler)->SetNumPreviousDirections(num_previous_directions); forest_valid = dynamic_cast*>(m_TrackingHandler)->IsForestValid(); } if (!forest_valid) { QMessageBox::information(nullptr, "Information", "Random forest is invalid. The forest signatue does not match the parameters of TrackingHandlerRandomForest."); return; } } } 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); } 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(m_Controls->m_SamplingDistanceBox->value()); tracker->SetUseStopVotes(m_Controls->m_StopVotesBox->isChecked()); tracker->SetOnlyForwardSamples(m_Controls->m_FrontalSamplesBox->isChecked()); tracker->SetAposterioriCurvCheck(false); tracker->SetMaxNumTracts(m_Controls->m_NumFibersBox->value()); tracker->SetNumberOfSamples(m_Controls->m_NumSamplesBox->value()); tracker->SetTrackingHandler(m_TrackingHandler); tracker->SetAngularThreshold(m_Controls->m_AngularThresholdBox->value()); tracker->SetMinTractLength(m_Controls->m_MinTractLengthBox->value()); tracker->SetUseOutputProbabilityMap(m_Controls->m_OutputProbMap->isChecked()); tracker->Update(); if (!tracker->GetUseOutputProbabilityMap()) { vtkSmartPointer fiberBundle = tracker->GetFiberPolyData(); if (!m_Controls->m_InteractiveBox->isChecked() && fiberBundle->GetNumberOfLines() == 0) { QMessageBox warnBox; warnBox.setWindowTitle("Warning"); warnBox.setText("No fiberbundle was generated!"); warnBox.setDetailedText("No fibers were generated using the chosen parameters. Typical reasons are:\n\n- Cutoff too high. Some images feature very low FA/GFA/peak size. Try to lower this parameter.\n- Angular threshold too strict. Try to increase this parameter.\n- A small step sizes also means many steps to go wrong. Especially in the case of probabilistic tractography. Try to adjust the angular threshold."); warnBox.setIcon(QMessageBox::Warning); warnBox.exec(); return; } mitk::FiberBundle::Pointer fib = mitk::FiberBundle::New(fiberBundle); fib->SetReferenceGeometry(dynamic_cast(m_InputImageNodes.at(0)->GetData())->GetGeometry()); if (m_Controls->m_ResampleFibersBox->isChecked()) fib->Compress(m_Controls->m_FiberErrorBox->value()); fib->ColorFibersByOrientation(); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(fib); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(fib); QString name("FiberBundle_"); name += m_InputImageNodes.at(0)->GetName().c_str(); name += "_Streamline"; node->SetName(name.toStdString()); GetDataStorage()->Add(node, m_InputImageNodes.at(0)); } } else { TrackerType::ItkDoubleImgType::Pointer outImg = tracker->GetOutputProbabilityMap(); mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(outImg.GetPointer()); img->SetVolume(outImg->GetBufferPointer()); if (m_Controls->m_InteractiveBox->isChecked()) { if (m_InteractiveNode.IsNull()) { m_InteractiveNode = mitk::DataNode::New(); QString name("Interactive"); m_InteractiveNode->SetName(name.toStdString()); GetDataStorage()->Add(m_InteractiveNode); } m_InteractiveNode->SetData(img); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); m_InteractiveNode->SetProperty("LookupTable", lut_prop); m_InteractiveNode->SetProperty("opacity", mitk::FloatProperty::New(0.5)); if (auto renderWindowPart = this->GetRenderWindowPart()) renderWindowPart->RequestUpdate(); } else { mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(img); QString name("ProbabilityMap_"); name += m_InputImageNodes.at(0)->GetName().c_str(); node->SetName(name.toStdString()); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType(mitk::LookupTable::JET_TRANSPARENT); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable(lut); node->SetProperty("LookupTable", lut_prop); node->SetProperty("opacity", mitk::FloatProperty::New(0.5)); GetDataStorage()->Add(node, m_InputImageNodes.at(0)); } } } 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 a6b065ddae..4304dcc8ac 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,1036 +1,1033 @@ QmitkStreamlineTrackingViewControls 0 0 413 1222 0 0 QmitkTemplate 0 0 3 3 Please Select Input Data Mask Image: Fibers that enter a region defined in this image will stop immediately. QComboBox::AdjustToMinimumContentsLength - Seed points are only placed inside the mask image. If no seed mask is selected, the whole image is seeded. QComboBox::AdjustToMinimumContentsLength - Stop Image: Tissue label image needed for gray matter seeding (WM=3, GM=1). Use e.g. MRtrix 5ttgen to generate such a label image. QComboBox::AdjustToMinimumContentsLength - Tractography is only performed inside the mask image. Fibers that leave the mask image are stopped. QComboBox::AdjustToMinimumContentsLength - FA/GFA Image: If an image is selected, the stopping criterion is not calculated from the input image but instead the selected image is used. QComboBox::AdjustToMinimumContentsLength - Seed Image: Tissue Image: Input Image. ODF, tensor, peak, and, in case of ML tractography, raw diffusion-weighted images are currently supported. <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Input Image. ODF, tensor and peak images are currently supported. Input Image: Tractography Forest: - Random forest for machine-learning based tractography. + Random forest for machine learning based tractography. QComboBox::AdjustToMinimumContentsLength - false Start Tractography 0 0 Interactive Tractography Number of seed points normally distributed around selected position. 1 9999999 50 Num. Seeds: true Dynamically pick seed location by click into image. Enable Interactive Tractography Seedpoints are normally distributed within a sphere centered at the selected position with the specified radius (in mm). 2 50.000000000000000 0.100000000000000 2.500000000000000 Radius: true When checked, parameter changes cause instant retracking while in interactive mode. Update on Parameter Change true 0 0 Parameters Mode: Toggle between deterministic and probabilistic tractography. Some modes might not be available for all types of tractography. Deterministic Probabilistic Seeds per Voxel: Number of seed points placed in each voxel. 1 9999999 Max. num. fibers: Tractography is stopped after the desired number of fibers is reached, even before all seed points are processed. -1 999999999 -1 Cutoff: Threshold on peak magnitude, FA, GFA, ... 5 1.000000000000000 0.100000000000000 0.100000000000000 ODF Cutoff: Additional threshold on the ODF magnitude. This is useful in case of CSD fODF tractography. 5 1.000000000000000 0.100000000000000 0.100000000000000 If you are using dODF images as input, it is advisable to sharpen the ODFs (min-max normalize and raise to the power of 4). This is not necessary for CSD fODFs, since they are naturally much sharper. Sharpen ODFs Qt::Horizontal QSizePolicy::Fixed 200 0 Advanced Parameters QFrame::NoFrame QFrame::Raised 0 0 0 0 Min. Tract Length: Step size (in voxels) 2 0.010000000000000 10.000000000000000 0.100000000000000 0.500000000000000 Angular Threshold: Flip directions: f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 1.000000000000000 f parameter of tensor tractography. f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). f: g: Shorter fibers are discarded. Minimum fiber length (in mm) 1 999.000000000000000 1.000000000000000 20.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 Internally flips progression directions. This might be necessary depending on the input data. x Internally flips progression directions. This might be necessary depending on the input data. y Internally flips progression directions. This might be necessary depending on the input data. z Step Size: Default: 90° * step_size -1 90 1 -1 f=1 + g=0 means FACT (depending on the chosen interpolation). f=0 and g=1 means TEND (disable interpolation for this mode!). 2 1.000000000000000 0.100000000000000 0.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 0 0 Requires tissue image. Enable Gray Matter Seeding false If false, nearest neighbor interpolation is used. Enable Trilinear Interpolation true 0 0 Neighborhood Sampling QFrame::NoFrame QFrame::Raised 0 0 0 0 Num. Samples: Number of neighborhood samples that are used to determine the next fiber progression direction. 50 Sampling Distance: Sampling distance (in voxels) 2 10.000000000000000 0.100000000000000 0.250000000000000 Only neighborhood samples in front of the current streamline position are considered. Use Only Frontal Samples false If checked, the majority of sampling points has to place a stop-vote for the streamline to terminate. If not checked, all sampling positions have to vote for a streamline termination. Use Stop-Votes false Qt::Vertical QSizePolicy::Expanding 20 220 0 0 Output and Postprocessing QFrame::NoFrame QFrame::Raised 0 0 0 0 Resample fibers using the specified error constraint. Compress Fibers true Qt::StrongFocus Lossy fiber compression. Recommended for large tractograms. Maximum error in mm. 3 10.000000000000000 0.010000000000000 0.100000000000000 Output probability map instead of tractogram. Output Probability Map false - m_InterpolationBox - frame_5 - m_OutputProbMap QmitkDataStorageComboBoxWithSelectNone QComboBox
QmitkDataStorageComboBoxWithSelectNone.h
 QmitkDataStorageComboBox QComboBox
QmitkDataStorageComboBox.h