diff --git a/Documentation/Doxygen/UserManual/MiniApps.dox b/Documentation/Doxygen/UserManual/MiniApps.dox new file mode 100644 index 0000000000..104522399f --- /dev/null +++ b/Documentation/Doxygen/UserManual/MiniApps.dox @@ -0,0 +1,70 @@ +/** +\page MiniAppExplainPage MITK MiniApps + +\section MiniAppDescription What are MiniApps + +MiniApps are small compilations of command line tools. Each of these tools is designed to fulfill one simple task, +e.g. resample an image or extract image statistics of a given region of interest (ROI). Several tools that relate to +a similiar topic or research area are grouped into one MiniApp. + +They are intended to provide command line access to a variety of features of MITK, thus facilitating batched processing of data. + +\section MiniAppUsage Usage + +The MiniApps are built in a self-describing way. When calling a MiniApp without any arguements it will list +all available sub-tools. When calling e.g. the DiffusionMiniApp it will look similarly to this: + +\code +$./MitkDiffusionMiniApps + +Please choose the mini app to execute: +(0) BatchedFolderRegistration +(1) CopyGeometry +(2) DicomFolderDump +(3) DiffusionIndices +(4) DwiDenoising +(5) ExportShImage +(6) ExtractImageStatistics +(7) FiberDirectionExtraction +(8) FiberProcessing +(9) FileFormatConverter +(10) GibbsTracking +(11) LocalDirectionalFiberPlausibility +(12) MultishellMethods +(13) NetworkCreation +(14) NetworkStatistics +(15) PeakExtraction +(16) PeaksAngularError +(17) QballReconstruction +(18) StreamlineTracking +(19) TensorDerivedMapsExtraction +(20) TensorReconstruction +Please select: +\endcode + +In order to select one of those tools simply append the displayed name to call, e.g. GibbsTracking +this will provide a listing of the parameters of that tool: +\code +$./MitkDiffusionMiniApps GibbsTracking + +[1.081] Start GibbsTracking .. + -i, --input, input image (tensor, Q-ball or FSL/MRTrix SH-coefficient image) + -p, --parameters, parameter file (.gtp) + -m, --mask, binary mask image (optional) + -s, --shConvention, sh coefficient convention (FSL, MRtrix) (optional), (default: FSL) + -o, --outFile, output fiber bundle (.fib) + -f, --noFlip, do not flip input image to match MITK coordinate convention (optional) +\endcode + +To execute the tool with parameters an exemplary call would look like this: + +\code +$./MitkDiffusionMiniApps GibbsTracking -i test.dti -p param.gtp -o /tmp/fiber.fib +\endcode + + +\section MiniAppAvailableList Available MiniApps + +\li \ref DiffusionMiniApps + +*/ diff --git a/Documentation/Doxygen/UserManual/UserManualPortal.dox b/Documentation/Doxygen/UserManual/UserManualPortal.dox index 995dddb941..94b5bf8e9d 100644 --- a/Documentation/Doxygen/UserManual/UserManualPortal.dox +++ b/Documentation/Doxygen/UserManual/UserManualPortal.dox @@ -1,20 +1,21 @@ /** \usersguidemainpage{UserManualPortal} MITK: User Manual To get an introduction to the usage of any MITK based application please read \ref MITKUserManualPage. It will give you an overview of most of the common questions, such as how to load or save data or navigate within it. This is a good starting point for first questions. Depending on what kind of work you intend do perform with MITK, certain applications are better suited to your needs than others. MITK offers a number of these Applications, each of which features a set of Plugins, which can solve certain tasks. To Learn more about MITK applications, please visit the \ref ApplicationsPage. For more specific information on how a plugin operates you can find the plugin documentation in \ref PluginListPage. The Plugin documentation usually explains the functionality in depth and should solve most problems you might encounter with the plugin. Depending on the application you are using you might have only some or all of the listed plugins available. Lastly, if your question is not answered here, please use our Mailinglist to let us know about your problem. Alternatively, you can contact us directly.

UserManualPortalTopics List of topics

-*/ \ No newline at end of file +*/ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h index 3a2c95f769..10e2d15222 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h @@ -1,230 +1,232 @@ /*=================================================================== 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_h_ #define _itk_TensorReconstructionWithEigenvalueCorrectionFilter_h_ #include "itkImageToImageFilter.h" #include #include #include #include namespace itk { template class TensorReconstructionWithEigenvalueCorrectionFilter : public ImageToImageFilter< itk::Image< TDiffusionPixelType, 3 >, itk::Image,3> > { public: typedef itk::VectorImage ImageType; typedef TensorReconstructionWithEigenvalueCorrectionFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< TDiffusionPixelType, 3>, Image< DiffusionTensor3D< TTensorPixelType >, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(TensorReconstructionWithEigenvalueCorrectionFilter, ImageToImageFilter) + typedef typename itk::Image,3> OutputType; typedef TDiffusionPixelType GradientPixelType; typedef DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef Image< TensorPixelType, 3 > TensorImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** An alternative typedef defining one (of the many) gradient images. * It will be assumed that the vectorImage has the same dimension as the * Reference image and a vector length parameter of \c n (number of * gradient directions) */ typedef VectorImage< GradientPixelType, 3 > GradientImagesType; typedef typename GradientImagesType::PixelType GradientVectorType; /** Holds each magnetic field gradient used to acquire one DWImage */ typedef vnl_vector_fixed< double, 3 > GradientDirectionType; /** Container to hold gradient directions of the 'n' DW measurements */ typedef VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; /** Another set method to add a gradient directions and its corresponding * image. The image here is a VectorImage. The user is expected to pass the * gradient directions in a container. The ith element of the container * corresponds to the gradient direction of the ith component image the * VectorImage. For the baseline image, a vector of all zeros * should be set. */ void SetGradientImage( GradientDirectionContainerType *, const GradientImagesType *image); + itkSetMacro( BValue, TTensorPixelType) /** Set the b0 threshold */ itkSetMacro( B0Threshold, double) /** Get the pseudeInverse that was calculated in the process of tensor estimation */ itkGetMacro(PseudoInverse, vnl_matrix) /** FIXME: added by Sebastian Wirkert due to compile error otherwise. */ itkGetMacro(CorrectedDiffusionVolumes, ImageType::Pointer) /** Outputs the design matrix that was calculated in the process of tensor estimation */ itkGetMacro(H, vnl_matrix) /** Outputs the normalized b-vector */ itkGetMacro(BVec, vnl_vector) /** Outputs the mask created by thresholding the B0 weighted volume*/ itkGetMacro(B0Mask, vnl_vector) /** Returns the dwi image volumes with smoothed voxels on positions were there were negative eigenvalues in the tensir*/ ImageType::Pointer GetGradientImagePointer() { return m_GradientImagePointer; } /** Get the mask image*/ itk::Image::Pointer GetMask() { return m_MaskImage; } protected: TensorReconstructionWithEigenvalueCorrectionFilter(); ~TensorReconstructionWithEigenvalueCorrectionFilter() {} void GenerateData(); typedef enum { GradientIsInASingleImage = 1, GradientIsInManyImages, Else } GradientImageTypeEnumeration; private: double CheckNeighbours(int x, int y, int z,int f, itk::Size<3> size, itk::Image ::Pointer mask,typename GradientImagesType::Pointer corrected_diffusion_temp); /** Calculates the attenuation for a voxel*/ void CalculateAttenuation(vnl_vector org_data, vnl_vector &atten,int nof,int numberb0); /** Correct the diffusion data set for voxels that contain negative eigenvalues in the tensor*/ void 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); /** Calculte a tensor iamge from a diffusion data set*/ void GenerateTensorImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,double what_mask, typename itk::Image< itk::DiffusionTensor3D, 3 >::Pointer tensorImg ); //void DeepCopyTensorImage(itk::Image< itk::DiffusionTensor3D, 3 >::Pointer tensorImg, itk::Image< itk::DiffusionTensor3D, 3 >::Pointer temp_tensorImg); //void DeepCopyDiffusionImage(itk::VectorImage::Pointer corrected_diffusion, itk::VectorImage::Pointer corrected_diffusion_temp,int nof); void TurnMask( itk::Size<3> size, itk::Image::Pointer mask, double previous_mask, double set_mask); double CheckNegatives ( itk::Size<3> size, itk::Image::Pointer mask, typename itk::Image< itk::DiffusionTensor3D, 3 >::Pointer tensorImg ); /** Gradient image was specified in a single image or in multiple images */ GradientImageTypeEnumeration m_GradientImageTypeEnumeration; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** b-value */ TTensorPixelType m_BValue; /** Number of baseline images */ unsigned int m_NumberOfBaselineImages; ImageType::Pointer m_CorrectedDiffusionVolumes; double m_B0Threshold; /** decodes what is to be done with every single voxel */ itk::Image::Pointer m_MaskImage; /** Pseudoinverse calculated in the process of calculating a tensor */ vnl_matrix m_PseudoInverse; /** design matrix */ vnl_matrix m_H; /** normalized b-vector */ vnl_vector m_BVec; /** m_B0Mask indicates whether a volume identified by an index is B0-weighted or not */ vnl_vector m_B0Mask; /** volumes of the diffusion data set */ typename GradientImagesType::Pointer m_GradientImagePointer; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkTensorReconstructionWithEigenvalueCorrectionFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx index 00f12b04e9..582b7ae4b9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx @@ -1,1098 +1,994 @@ /*=================================================================== 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_TensorReconstructionWithEigenvalueCorrectioFiltern_txx_ +#define _itk_TensorReconstructionWithEigenvalueCorrectionFilter_txx_ #endif #include "itkImageRegionConstIterator.h" #include #include "itkImageFileWriter.h" #include "itkImage.h" #include "itkImageRegionIterator.h" - - namespace itk { +template +TensorReconstructionWithEigenvalueCorrectionFilter +::TensorReconstructionWithEigenvalueCorrectionFilter() +{ + m_B0Threshold = 50.0; +} - template - TensorReconstructionWithEigenvalueCorrectionFilter - ::TensorReconstructionWithEigenvalueCorrectionFilter() - { - m_B0Threshold = 50.0; - } - - template - void - TensorReconstructionWithEigenvalueCorrectionFilter - ::GenerateData () - { +template +void +TensorReconstructionWithEigenvalueCorrectionFilter +::GenerateData () +{ - m_GradientImagePointer = static_cast< GradientImagesType * >( - this->ProcessObject::GetInput(0) ); + m_GradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); - typename GradientImagesType::SizeType size = m_GradientImagePointer->GetLargestPossibleRegion().GetSize(); + typename GradientImagesType::SizeType size = m_GradientImagePointer->GetLargestPossibleRegion().GetSize(); - // number of volumes - int nof = m_GradientDirectionContainer->Size(); + // 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); + // due to roundings, the values are not always exactly zero + if(vec[0]<0.0001 && vec[1]<0.0001 && vec[2]<0.0001 && vec[0]>-0.0001&& vec[1]>-0.0001 && vec[2]>-0.0001) { - vnl_vector_fixed vec = m_GradientDirectionContainer->ElementAt(i); - - // due to roundings, the values are not always exactly zero - if(vec[0]<0.0001 && vec[1]<0.0001 && vec[2]<0.0001 && vec[0]>-0.0001&& vec[1]>-0.0001 && vec[2]>-0.0001) - { - numberb0++; - } + numberb0++; } + } - // Matrix to store all diffusion encoding gradients - vnl_matrix directions(nof-numberb0,3); + // Matrix to store all diffusion encoding gradients + vnl_matrix directions(nof-numberb0,3); - m_B0Mask.set_size(nof); + m_B0Mask.set_size(nof); - int cnt=0; - for(int i=0; i vec = m_GradientDirectionContainer->ElementAt(i); + int cnt=0; + for(int i=0; i vec = m_GradientDirectionContainer->ElementAt(i); - if(vec[0]<0.0001 && vec[1]<0.0001 && vec[2]<0.0001 && vec[0]>-0.0001&& vec[1]>-0.0001 && vec[2]>-0.0001) - { - // 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; + if(vec[0]<0.0001 && vec[1]<0.0001 && vec[2]<0.0001 && vec[0]>-0.0001&& vec[1]>-0.0001 && vec[2]>-0.0001) + { + // 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]; + // 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++; - } + cnt++; } + } - // looking for maximal norm among gradients. - // The norm is calculated with use of spectral radius theorem- based on determination of eigenvalue. + // 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_matrix dirsTimesDirsTrans = directions*directions.transpose(); - vnl_vector< double> diagonal(nof-numberb0); - vnl_vector< double> b_vec(nof-numberb0); - vnl_vector< double> temporary(3); + 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); + // Calculation of so called design matrix that is used to obtain expected signal. - //H is matrix that containes covariances for directions. It is stored twice because its original value is needed later - // while H is changed + vnl_matrix H(nof-numberb0, 6); + vnl_matrix H_org(nof-numberb0, 6); + vnl_vector pre_tensor(9); - int etbt[6] = { 0, 4, 8, 1, 5, 2 };// tensor order + //H is matrix that contains covariances for directions. It is stored twice because its original value is needed later + // while H is changed - for (int i = 0; i < nof-numberb0; i++) + 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++) { - 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++) + temporary[j] = -directions[i][j]; + } + for (int j = 0; j < 3; j++) + { + for (int k = 0; k < 3; k++) { - H[i][3 + j] *= 2.0; + 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; + H_org=H; - // calculation of inverse matrix by means of pseudoinverse + // 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(); + 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(); + 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. + // 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; - for(unsigned int x=0;x ix = {{x,y,z}}; - + double mean_b=0.0; + itk::Index<3> ix = {{x,y,z}}; - GradientVectorType pixel = m_GradientImagePointer->GetPixel(ix); - for (int i=0;i m_B0Threshold) - { - mask->SetPixel(ix, 1); - mask_cnt++; - } - else + GradientVectorType pixel = m_GradientImagePointer->GetPixel(ix); + for (int i=0;iSetPixel(ix, 0); + mean_b = mean_b + pixel[i]; } - } - } - } - - //14.10.2013 - /* - - //create a copy of the original image- it is then used in pre-processing methods - - m_CorrectedDiffusionVolumes = ImageType::New(); - m_CorrectedDiffusionVolumes->SetRegions(size); - m_CorrectedDiffusionVolumes->SetSpacing(m_GradientImagePointer->GetSpacing()); - m_CorrectedDiffusionVolumes->SetOrigin(m_GradientImagePointer->GetOrigin()); - m_CorrectedDiffusionVolumes->SetVectorLength(nof); - m_CorrectedDiffusionVolumes->Allocate(); - - - - for ( int x=0;x m_B0Threshold) { - itk::Index<3> ix = {{x,y,z}}; - - GradientVectorType p = m_GradientImagePointer->GetPixel(ix); - - m_CorrectedDiffusionVolumes->SetPixel(ix,p); - + mask->SetPixel(ix, 1); + mask_cnt++; } + else + mask->SetPixel(ix, 0); } } + } - */ - - - - //Sometimes the gradient voxels may contain negative values ( even if B0 voxel is > = 50 ). This must be corrected by smoothing DWI - //Smoothing is done by aproximation of negative voxel value by its correct ( positive) 27-th neighborhood. - - - // typename TensorImageType::Pointer someimg; - // someimg = TensorImageType::New(); - - - - double mask_val=0.0; - vnl_vector org_vec(nof); + double mask_val=0.0; + vnl_vector org_vec(nof); - int counter_corrected =0; + int counter_corrected =0; - for (unsigned int x=0;x ix = {x,y,z}; + itk::Index<3> ix = {x,y,z}; - mask_val = mask->GetPixel(ix); + mask_val = mask->GetPixel(ix); - //14.10.2013 - //GradientVectorType pixel2 = m_CorrectedDiffusionVolumes->GetPixel(ix); + GradientVectorType pixel2 = m_GradientImagePointer->GetPixel(ix); - GradientVectorType pixel2 = m_GradientImagePointer->GetPixel(ix); - - for (int i=0;i0) + if(mask_val >0) + { + for( int f=0;fSetPixel(ix, pixel2); - m_GradientImagePointer->SetPixel(ix, pixel2); + for (int i=0;iSetPixel(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 - 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. + for (int i=0;i pixel_max(nof-numberb0); - vnl_vector< double> pixel_min(nof-numberb0); + m_PseudoInverse = pseudoInverse; + m_H = H_org; + m_BVec=b_vec; - // to high and to low attenuation is calculated with use of highest allowed =5 and lowest allowed =0.01 diffusion coefficient + // in preprocessing we are dealing with 3 types of voxels: voxels excluded = 0 in mask, voxels correct =1 and voxels under correction + //= 2. During pre processing most of voxels should be switched from 2 to 1. - for (int i=0;i outputIterator(tensorImg, tensorImg->GetLargestPossibleRegion()); + outputIterator.GoToBegin(); + while(!outputIterator.IsAtEnd()) + { + TensorPixelType tens = outputIterator.Get(); - // Generation of final pre-processed tensor image that might be used as an input for FWE method - //14.10.2013 - //GenerateTensorImage(nof,numberb0,size,m_CorrectedDiffusionVolumes,mask,what_mask,tensorImg); - GenerateTensorImage(nof,numberb0,size,m_GradientImagePointer,mask,what_mask,tensorImg); + tens/= 1000.0; - m_MaskImage = mask; + outputIterator.Set(tens); + ++outputIterator; + } - this->SetNthOutput(0, tensorImg); + this->SetNthOutput(0, tensorImg); +} - } +template +void +TensorReconstructionWithEigenvalueCorrectionFilter +::SetGradientImage( GradientDirectionContainerType *gradientDirection, + const GradientImagesType *gradientImage ) +{ - 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."); + } - if( m_GradientImageTypeEnumeration == GradientIsInManyImages ) - { - itkExceptionMacro( << "Cannot call both methods:" - << "AddGradientImage and SetGradientImage. Please call only one of them."); - } - - this->m_GradientDirectionContainer = gradientDirection; - + this->m_GradientDirectionContainer = gradientDirection; - unsigned int numImages = gradientDirection->Size(); - this->m_NumberOfBaselineImages = 0; - this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; + unsigned int numImages = gradientDirection->Size(); + this->m_NumberOfBaselineImages = 0; - // 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->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages; - this->ProcessObject::SetNthInput( 0, - const_cast< GradientImagesType* >(gradientImage) ); - m_GradientImageTypeEnumeration = GradientIsInASingleImage; + // 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."); } - 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. + 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; + // 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 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 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; + double tempsum=0; + double temp_number=0; - for(int i=back_x; i<=forth_x; i++) + for(int i=back_x; i<=forth_x; i++) + { + for (int j=back_y; j<=forth_y; j++) { - for (int j=back_y; j<=forth_y; j++) + for (int k=back_z; k<=forth_z; k++) { - for (int k=back_z; k<=forth_z; k++) - { - itk::Index<3> ix = {i,j,k}; + itk::Index<3> ix = {i,j,k}; - GradientVectorType p = corrected_diffusion_temp->GetPixel(ix); + GradientVectorType p = corrected_diffusion_temp->GetPixel(ix); - if (p[f] > 0.0 )// taking only positive values and counting them + if (p[f] > 0.0 )// taking only positive values and counting them + { + if(!(i==x && j==y && k== z)) { - if(!(i==x && j==y && k== z)) - { - tempsum=tempsum+p[f]; - temp_number++; - } - + tempsum=tempsum+p[f]; + temp_number++; } - - } } } + } - //getting back to the original position of the voxel + //getting back to the original position of the voxel - itk::Index<3> ix = {x,y,z}; + itk::Index<3> ix = {x,y,z}; - if (temp_number <= 0.0) - { - tempsum=0; - mask->SetPixel(ix,0); - } - else - { - tempsum=tempsum/temp_number; + if (temp_number <= 0.0) + { + tempsum=0; + mask->SetPixel(ix,0); + } + else + { + tempsum=tempsum/temp_number; - } + } - return tempsum;// smoothed value of voxel - } + 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; +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) { - if(m_B0Mask[i]>0) - { - double o_d=org_data[i]; - mean_b=mean_b+org_data[i]; - } - + double o_d=org_data[i]; + 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 ) - { +template +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. - // 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; + double pixel=0; + itk::DiffusionTensor3D ten; + vnl_matrix temp_tensor(3,3); + vnl_vector eigen_vals(3); + vnl_vector tensor (6); - // declaration of important structures and variables - double badvoxels=0; - double pixel=0; - itk::DiffusionTensor3D ten; - vnl_matrix temp_tensor(3,3); - vnl_vector eigen_vals(3); - vnl_vector tensor (6); + // for every pixel from the image + for (unsigned int x=0;x ix = {x,y,z}; + pixel = mask->GetPixel(ix); + // but only if previously marked as bad one-negative eigen value - itk::Index<3> ix = {x,y,z}; - pixel = mask->GetPixel(ix); + if(pixel > 1) + { - // but only if previously marked as bad one-negative eigen value + ten = tensorImg->GetPixel(ix); - if(pixel > 1) - { + // changing order from tensor structure in MITK into vnl structure 3x3 symmetric tensor matrix - ten = tensorImg->GetPixel(ix); + tensor[0] = ten(0,0); + tensor[3] = ten(0,1); + tensor[5] = ten(0,2); + tensor[1] = ten(1,1); + tensor[4] = ten(1,2); + tensor[2] = ten(2,2); - // changing order from tensor structure in MITK into vnl structure 3x3 symmetric tensor matrix + temp_tensor[0][0]= tensor[0]; temp_tensor[1][0]= tensor[3]; temp_tensor[2][0]= tensor[5]; + temp_tensor[0][1]= tensor[3]; temp_tensor[1][1]= tensor[1]; temp_tensor[2][1]= tensor[4]; + temp_tensor[0][2]= tensor[5]; temp_tensor[1][2]= tensor[4]; temp_tensor[2][2]= tensor[2]; - tensor[0] = ten(0,0); - tensor[3] = ten(0,1); - tensor[5] = ten(0,2); - tensor[1] = ten(1,1); - tensor[4] = ten(1,2); - tensor[2] = ten(2,2); + //checking negativity of tensor eigenvalues - temp_tensor[0][0]= tensor[0]; temp_tensor[1][0]= tensor[3]; temp_tensor[2][0]= tensor[5]; - temp_tensor[0][1]= tensor[3]; temp_tensor[1][1]= tensor[1]; temp_tensor[2][1]= tensor[4]; - temp_tensor[0][2]= tensor[5]; temp_tensor[1][2]= tensor[4]; temp_tensor[2][2]= tensor[2]; + vnl_symmetric_eigensystem eigen_tensor(temp_tensor); - //checking negativity of tensor eigenvalues - vnl_symmetric_eigensystem eigen_tensor(temp_tensor); + eigen_vals[0]=eigen_tensor.get_eigenvalue(0); + eigen_vals[1]=eigen_tensor.get_eigenvalue(1); + eigen_vals[2]=eigen_tensor.get_eigenvalue(2); + //comparison to 0.01 instead of 0 was proposed by O.Pasternak - eigen_vals[0]=eigen_tensor.get_eigenvalue(0); - eigen_vals[1]=eigen_tensor.get_eigenvalue(1); - eigen_vals[2]=eigen_tensor.get_eigenvalue(2); + if( eigen_vals[0]>0.01 && eigen_vals[1]>0.01 && eigen_vals[2]>0.01) + { + mask->SetPixel(ix,1); - //comparison to 0.01 instead of 0 was proposed by O.Pasternak + } + else + { + badvoxels++; + } - if( eigen_vals[0]>0.01 && eigen_vals[1]>0.01 && eigen_vals[2]>0.01) - { - mask->SetPixel(ix,1); + } - } - else - { - badvoxels++; - } + } + } + } - } + double ret = badvoxels; - } - } - } + return ret; - double ret = badvoxels; +} - return ret; - } +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 + vnl_vector org_data(nof); + vnl_vector atten(nof-numberb0); + double cnt_atten=0; - 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) + for (unsigned int z=0;z org_data(nof); - vnl_vector atten(nof-numberb0); - double cnt_atten=0; - - for (unsigned int z=0;z ix = {x, y, z}; - for (unsigned int y=0;y ix = {x, y, z}; + if(mask->GetPixel(ix) > 1.0) + { + GradientVectorType pt = corrected_diffusion->GetPixel(ix); - if(mask->GetPixel(ix) > 1.0) + for (int i=0;iGetPixel(ix); - - for (int i=0;i0) { - if(m_B0Mask[i]>0) - { - mean_b=mean_b+org_data[i]; - } + 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; + 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 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 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; + double tempsum=0; + double temp_number=0; - for(unsigned int i=back_x; i<=forth_x; i++) + for(unsigned int i=back_x; i<=forth_x; i++) + { + for (unsigned int j=back_y; j<=forth_y; j++) { - for (unsigned int j=back_y; j<=forth_y; j++) + for (unsigned int k=back_z; k<=forth_z; k++) { - for (unsigned int k=back_z; k<=forth_z; k++) - { - itk::Index<3> ix = {i,j,k}; + itk::Index<3> ix = {i,j,k}; - GradientVectorType p = corrected_diffusion->GetPixel(ix); + GradientVectorType p = corrected_diffusion->GetPixel(ix); - //double test= p[f]; + //double test= p[f]; - if (p[f] > 0.0 )// taking only positive values and counting them + if (p[f] > 0.0 )// taking only positive values and counting them + { + if(!(i==x && j==y && k== z)) { - if(!(i==x && j==y && k== z)) - { - tempsum=tempsum+p[f]; - temp_number++; - } - + tempsum=tempsum+p[f]; + temp_number++; } - } + + } } + } - //getting back to the original position of the voxel + //getting back to the original position of the voxel - itk::Index<3> ix = {x,y,z}; + itk::Index<3> ix = {x,y,z}; - if (temp_number <= 0.0) - { - tempsum=0; - mask->SetPixel(ix,0); - } - else - { - tempsum=tempsum/temp_number; + if (temp_number <= 0.0) + { + tempsum=0; + mask->SetPixel(ix,0); + } + else + { + tempsum=tempsum/temp_number; - } + } - org_data[f] = tempsum; + org_data[f] = tempsum; } cnt_atten++; - } - //smoothing B0 - if(m_B0Mask[f]==1) - { + } + //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; + 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 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 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; + double tempsum=0; + double temp_number=0; - for(unsigned int i=back_x; i<=forth_x; i++) + 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 (unsigned int j=back_y; j<=forth_y; j++) - { - for (unsigned int k=back_z; k<=forth_z; k++) - { - itk::Index<3> ix = {i,j,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++; - } - - } + 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 = {x,y,z}; - if (temp_number <= 0.0) - { - tempsum=0; - mask->SetPixel(ix,0); } - else - { - tempsum=tempsum/temp_number; + } + } - } + //getting back to the original position of the voxel - org_data[f] = tempsum; + itk::Index<3> ix = {x,y,z}; + if (temp_number <= 0.0) + { + tempsum=0; + 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); - corrected_diffusion->SetPixel(ix, pt); + for (int i=0;iSetPixel(ix, pt); + + } + else + { + GradientVectorType pt = corrected_diffusion->GetPixel(ix); + 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 what_mask, 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); +template +void +TensorReconstructionWithEigenvalueCorrectionFilter +::GenerateTensorImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,double what_mask, 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); - itk::Index<3> ix; - vnl_vector org_data(nof); - vnl_vector atten(nof-numberb0); - vnl_vector tensor(6); - itk::DiffusionTensor3D ten; - double mask_val=0; + itk::Index<3> ix; + itk::Index<3> idx; + idx.Fill(5); - for (unsigned int x=0;x org_data(nof); + vnl_vector atten(nof-numberb0); + vnl_vector tensor(6); + itk::DiffusionTensor3D ten; + double mask_val=0; - for (unsigned int y=0;yGetPixel(ix); - //Tensors are calculated only for voxels above theshold for B0 image. + ix[0] = x; ix[1] = y; ix[2] = z; - if( mask_val > 0.0 ) - { + mask_val= mask->GetPixel(ix); - // calculation of attenuation with use of gradient image and and mean B0 image - GradientVectorType pt = corrected_diffusion->GetPixel(ix); + //Tensors are calculated only for voxels above theshold for B0 image. - for (int i=0;i 0.0 ) + { - double mean_b=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]; - } + 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); + // Calculation of tensor with use of previously calculated inverse of design matrix and attenuation + tensor = m_PseudoInverse*atten; + ten(0,0) = tensor[0]; + ten(0,1) = tensor[3]; + ten(0,2) = tensor[5]; + ten(1,1) = tensor[1]; + ten(1,2) = tensor[4]; + ten(2,2) = tensor[2]; + tensorImg->SetPixel(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; - tensorImg->SetPixel(ix, ten); - } - } } - } - - + // for voxels with mask value 0 - tensor is simply 0 ( outside brain value) + else if (mask_val < 1.0) + { - }// end of Generate Tensor + ten(0,0) = 0; + ten(0,1) = 0; + ten(0,2) = 0; + ten(1,1) = 0; + ten(1,2) = 0; + ten(2,2) = 0; + tensorImg->SetPixel(ix, ten); + } + if (ix == idx) + { + for (int ll = 0; ll < 6 ; ll++) + std::cout << tensor[ll] << "," << std::endl; + } - 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; - for(unsigned int x=0;xGetPixel(ix); - if(temp_mask_value>previous_mask) - { - mask->SetPixel(ix,set_mask); - } - } } } } - /* +}// end of Generate Tensor - template - void - TensorReconstructionWithEigenvalueCorrectionFilter - ::DeepCopyDiffusionImage(itk::VectorImage::Pointer corrected_diffusion, itk::VectorImage::Pointer corrected_diffusion_temp,int nof) - { - corrected_diffusion_temp->SetSpacing(corrected_diffusion->GetSpacing()); - corrected_diffusion_temp->SetOrigin(corrected_diffusion->GetOrigin()); - corrected_diffusion_temp->SetVectorLength(nof); - corrected_diffusion_temp->SetRegions(corrected_diffusion->GetLargestPossibleRegion()); - corrected_diffusion_temp->Allocate(); - itk::ImageRegionConstIterator inputIterator(corrected_diffusion, corrected_diffusion->GetLargestPossibleRegion()); - itk::ImageRegionIterator outputIterator(corrected_diffusion_temp, corrected_diffusion_temp->GetLargestPossibleRegion()); +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. - inputIterator.GoToBegin(); - outputIterator.GoToBegin(); + itk::Index<3> ix; + double temp_mask_value=0; - while(!inputIterator.IsAtEnd()) - { - outputIterator.Set(inputIterator.Get()); - ++inputIterator; - ++outputIterator; - } - } - - template - void - TensorReconstructionWithEigenvalueCorrectionFilter - ::DeepCopyTensorImage(itk::Image< itk::DiffusionTensor3D, 3 >::Pointer tensorImg, itk::Image< itk::DiffusionTensor3D, 3 >::Pointer temp_tensorImg) + for(unsigned int x=0;xSetSpacing(tensorImg->GetSpacing()); - temp_tensorImg->SetOrigin(tensorImg->GetOrigin()); - temp_tensorImg->SetRegions(tensorImg->GetLargestPossibleRegion()); - temp_tensorImg->Allocate(); - - itk::ImageRegionConstIterator inputIterator(tensorImg, tensorImg->GetLargestPossibleRegion()); - itk::ImageRegionIterator outputIterator(temp_tensorImg, temp_tensorImg->GetLargestPossibleRegion()); - - inputIterator.GoToBegin(); - outputIterator.GoToBegin(); - - while(!inputIterator.IsAtEnd()) + for(unsigned int y=0;yGetPixel(ix); + + if(temp_mask_value>previous_mask) + { + mask->SetPixel(ix,set_mask); + } } + } } - - - */ - - +} } // end of namespace diff --git a/Modules/DiffusionImaging/MiniApps/Documentation/DiffusionMiniApps.dox b/Modules/DiffusionImaging/MiniApps/Documentation/DiffusionMiniApps.dox new file mode 100644 index 0000000000..36ba822295 --- /dev/null +++ b/Modules/DiffusionImaging/MiniApps/Documentation/DiffusionMiniApps.dox @@ -0,0 +1,180 @@ +/** +\page DiffusionMiniApps MITK Diffusion MiniApps + +\tableofcontents + +The respective MiniApp is called MitkDiffusionMiniApp and is shipped with the current MITK Diffusion installer. +This page intends to provide an overview of all tools that are included in the DiffusionMiniApp. Also it relates them +to the respective Plugin in the MITK Diffusion application (if one exists). +For a detailed list of parameters call the according tool without any arguments (see \ref MiniAppExplainPage for details on this) or refer its Plugin equivalent. + +\section Preprocessing Preprocessing Tools + + +\subsection BatchedFolderRegistratuion Batched Folder Registration + +Allows to register a series of images (of different modalities, including diffusion weighted) to one reference image. It allows to register derived resources (e.g. a segmentation on +a T2 image) using the transformation of the original (T2) image. + +For the following examples assume a folder containing a longitudinal study with T1,T2, DWI images and segmentations (ROI) : + +\code +Patien01_2010-1.dwi +Patien01_2010-1_T1.nrrd +Patien01_2010-1_T2.nrrd +Patien01_2010-1_ROI.nrrd +Patien01_2010-2.dwi +Patien01_2010-2_T1.nrrd +Patien01_2010-2_T2.nrrd +Patien01_2010-2_ROI.nrrd +Patien01_2010-3.dwi +Patien01_2010-3_T1.nrrd +Patien01_2010-3_T2.nrrd +Patien01_2010-3_ROI.nrrd +Patien01_2010-4.dwi +Patien01_2010-4_T1.nrrd +Patien01_2010-4_T2.nrrd +Patien01_2010-4_ROI.nrrd +\endcode + +All T2 and DWI images are to be co-registered to the first T2 image, this can be achieved by the following two calls: + +\code + $./MitkDiffusionMiniApps BatchedFolderRegistration -i /home/inputFolder/ -o /home/outputFolder/ -f Patien01_2010-1_T2.nrrd -m T2.nrrd + $./MitkDiffusionMiniApps BatchedFolderRegistration -i /home/inputFolder/ -o /home/outputFolder/ -f Patien01_2010-1_T2.nrrd -m .dwi +\endcode + +The segmentations where performed on the T1 image and are therefore related to the image space of the respective T1 image, +so they can be bound to these images by marking them as derived resources. To register them both you would call + +\code +$./bin/MitkDiffusionMiniApps BatchedFolderRegistration -i /home/inputFolder/ -o /home/outputFolder/ -f Patien01_2010-1_T2.nrrd -m _T1.nrrd -d _ROI.nrrd -b +\endcode + +\note the suffixes of '_T1.nrrd' and '_ROI.nrrd' must have the same length! + +The parameter -b designates the derived resource as binary such that a nearest neighbor interpolation is used. + +All images (execpt for DWI files) are resample to the reference image, to resample to a specific spacing append the desired +spacing like this (e.g. 1 x 1 x 2 mm) + +\code +$./MitkDiffusionMiniApps BatchedFolderRegistration -i /home/inputFolder/ -o /home/outputFolder/ -f Patien01_2010-1_T2.nrrd -m .dwi -r 1,1,2 +\endcode + +\note Registration methods assume that both images occupy roughly the same space. It may happend that this is not the case, +and therefore registration fails. In this case you can try the -c option which uses the same origin for both images. + + +\subsection CopyGeometry Copy Geometry + +Copies the geometry (origin) of the source image to the target image. + +\subsection DicomLoad Dicom Loader + +Dicom Tools allow to parse dicom folders and export NRRD or DWI files, using standard naming. + +TODO enhance Docu, when MiniApp is ready .. + +\subsection TensorRecon Tensor Reconstruction +See \ref QmitkDiffusionImagingUserManualTensorReconstruction for the GUI equivalent of this tool. + +Takes a .dwi, .fsl/.fslgz file as input and saves the computed reconstructed tensor to the specified file. +It also allows for a threshold to be set, to exclude low b values from the reconstruction process. + +\code +./MitkDiffusionMiniApps TensorReconstruction -i /home/user/sample.dwi -o /home/user/tensors.dti -t 50 +\endcode + +\subsection QballRecon Qball Reconstruction + +See \ref QmitkDiffusionImagingUserManualQBallReconstruction for the GUI equivalent of this tool. + +\code +./MitkDiffusionMiniApps QballReconstruction -i /home/user/sample.dwi -o /home/user/tensors.qbi -t 50 -r .006 -shc /home/user/coeffs.csv +\endcode + +\subsection PeakExtraction Peak Extraction + + +\subsection PeakAngularErr Peak Angular Error + +\section DiffusionMeasures Diffusion Related Measures + +\subsection DiffusionIndices Diffusion Indices + +See \ref QmitkDiffusionImagingUserManualQuantification for the GUI equivalent of this tool. + +Computes a selected tensor derived indices (fa, gfa, ra, ad, rd, ca, l2, l3, md) given a +Tensor, Q-ball or FSL/MRTrix SH-coefficient image. E.g. to compute the fraction anisotropy call + +\code +./MitkDiffusionMiniApps DiffusionIndices -i /home/user/input.dti -idx fa -o /home/user/fa_image.nrrd +\endcode + +\subsection AllDiffusionIndices Tensor Derived Maps Extraction + +Similar to \ref DiffusionIndices . But computes all of the following indices FA, RA, MD, CA, RD, AD at once. +Also the input is a regular .dwi file, the tensor reconstruction is done implicitly (using a b0 threshold of 50). + + +\section FibTracking Fiber Tracking and Processing Methods + +\subsection FibDirection Fiber Direction Extraction + +TODO Peter ? + +\subsection Streamline Streamline Tracking + +See \ref org_mitk_views_streamlinetracking for the GUI equivalent of this tool. + +Performs a streamline tracking on a tensor image. + +\subsection GibbsTracking Gibbs Fiber Tracking + +See \ref org_mitk_views_gibbstracking for the GUI equivalent of this tool. + +Performs a Gibbs tracking on a tensor image. + +\subsection FibProcessing Fiber Processing + +Post-process a fiber bundle. Provides the possibility to + +\li remove short/long fiber tracks +\li combine fiber bundles +\li resample a fiber bundle +\li scale bundle in each direction independently + + +\subsection FibFoxProcessing Fiber Fox Processing + +See \ref QmitkFiberfoxViewUserManualSignalGeneration for the GUI equivalent of this tool. + +Generates a signal from a fiber bundle provided a reference DWI and a parameter file. The parameter file can be generated +using the Fiberfox plugin (sub-tab) Signal Generation. + +\subsection FormatConv File Format Converter + +Determines the data type and converts the input file (if possible) to .NRRD (regular image), +.DWI (diffusion image) or .FIB (fiber bundle). + +\subsection MultiShell Multishell Methods + +Computes several fits on an images (Kurtosis,Bi-Exponential, ADC). + +These fits are part of the Preprocessing Plugin \ref QmitkDiffusionImagingUserManualPreprocessing . + +\section NetworkTools Connectomics + +\subsection NetworkCreation Network Creation + +See \ref org_mitk_views_connectomicsdata for the GUI equivalent of this tool. + +Creates a network based on a brain parcellation and a fiber image. + +\subsection NetworkStatistics Network Statistics + +See \ref org_mitk_views_connectomicsstatistics for the GUI equivalent of this tool. + +Calculates several network statistics for a given connectome. +*/ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui index ef12cfca5d..cd7b9890cc 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui @@ -1,562 +1,562 @@ QmitkTensorReconstructionViewControls 0 0 380 1019 0 0 true QmitkTensorReconstructionViewControls Please Select Input Data: <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true Raw DWI/DTI: Tensor reconstruction and colormap visualization Advanced Settings false QFrame::StyledPanel QFrame::Raised 9 9 9 9 QFrame::NoFrame QFrame::Raised 0 0 0 0 B0 Threshold false 10000 Only influences WLS reconstruction Ignore voxels with negative eigenvalues 0 - Weighted Linear Least Squares + ITK Linear Least Squares With correction for negative eigenvalues false Select raw DWI! Start Reconstruction Estimate Diffusion Weighted Image from Tensors QFrame::NoFrame QFrame::Raised QFormLayout::AllNonFixedFieldsGrow 6 6 0 0 0 0 how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" B-Value false #Gradient Directions 3 12 42 92 162 252 362 492 642 812 1002 10000 100 1000 false Estimates the original diffusion weighted image based on a reconstructed tensor image. Estimate DWI based on Tensor Image Estimate Q-Ball Image from Tensors false Calculate ODF value as tensor value in the according direction Start QBI Estimation Estimate Residuals false false false percentages of error 1 Per volume 200 300 Per slice outliers per slice QFrame::NoFrame QFrame::Raised 0 0 0 0 300 400 QFrame::NoFrame QFrame::Raised 0 0 0 0 20 255 Volume: .., Slice:.. false Calculate the residual from a dti and a dwi image Residual Image Calculation QmitkResidualAnalysisWidget QWidget
QmitkResidualAnalysisWidget.h
1 QmitkResidualViewWidget QGraphicsView
QmitkResidualViewWidget.h