diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h index 3f2c6c5fb0..d89cd634ff 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h @@ -1,270 +1,203 @@ /*=================================================================== 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. ===================================================================*/ -/*========================================================================= - -Program: Tensor ToolKit - TTK -Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.h $ -Language: C++ -Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ -Version: $Revision: 68 $ - -Copyright (c) INRIA 2010. All rights reserved. -See LICENSE.txt for details. - -This software is distributed WITHOUT ANY WARRANTY; without even -the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the above copyright notices for more information. - -=========================================================================*/ #ifndef _itk_TensorReconstructionWithEigenvalueCorrectionFilter_h_ #define _itk_TensorReconstructionWithEigenvalueCorrectionFilter_h_ #include "itkImageToImageFilter.h" #include #include #include #include -#include - - - - -typedef itk::VectorImage ImageType; - 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. */ - itkNewMacro(Self); + itkNewMacro(Self) /** Runtime information support. */ - itkTypeMacro(TensorReconstructionWithEigenvalueCorrectionFilter, ImageToImageFilter); + itkTypeMacro(TensorReconstructionWithEigenvalueCorrectionFilter, ImageToImageFilter) + - typedef TDiffusionPixelType ReferencePixelType; typedef TDiffusionPixelType GradientPixelType; typedef DiffusionTensor3D< TTensorPixelType > TensorPixelType; - - /** Reference image data, This image is aquired in the absence - * of a diffusion sensitizing field gradient */ - typedef typename Superclass::InputImageType ReferenceImageType; typedef Image< TensorPixelType, 3 > TensorImageType; - typedef TensorImageType OutputImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - /** Typedef defining one (of the many) gradient images. */ - typedef Image< GradientPixelType, 3 > GradientImageType; /** 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 the tensor basis coefficients G_k - typedef vnl_matrix_fixed< double, 6, 6 > TensorBasisMatrixType; - - typedef vnl_matrix< double > CoefficientMatrixType; - */ /** 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); - /** Set method to set the reference image. */ - void SetReferenceImage( ReferenceImageType *referenceImage ) - { - if( m_GradientImageTypeEnumeration == GradientIsInASingleImage) - { - itkExceptionMacro( << "Cannot call both methods:" - << "AddGradientImage and SetGradientImage. Please call only one of them."); - } - this->ProcessObject::SetNthInput( 0, referenceImage ); - m_GradientImageTypeEnumeration = GradientIsInManyImages; - } - /** Get reference image */ - virtual ReferenceImageType * GetReferenceImage() - { return ( static_cast< ReferenceImageType *>(this->ProcessObject::GetInput(0)) ); } + itkSetMacro( BValue, TTensorPixelType) + itkSetMacro( B0Threshold, float) - /** Return the gradient direction. idx is 0 based */ - virtual GradientDirectionType GetGradientDirection( unsigned int idx) const - { - if( idx >= m_NumberOfGradientDirections ) - { - itkExceptionMacro( << "Gradient direction " << idx << "does not exist" ); - } - return m_GradientDirectionContainer->ElementAt( idx+1 ); - } - - itkSetMacro( BValue, TTensorPixelType); - itkSetMacro( B0Threshold, float); - itkSetMacro (Flagstatus, int); + itkGetMacro(PseudoInverse, vnl_matrix) + itkGetMacro(H, vnl_matrix) + itkGetMacro(BVec, vnl_vector) + itkGetMacro(B0Mask, vnl_vector) - itkGetMacro(PseudoInverse, vnl_matrix); - itkGetMacro(H, vnl_matrix); - itkGetMacro(BVec, vnl_vector); - itkGetMacro(B0Mask, vnl_vector); - itkGetMacro(Voxdim, vnl_vector); - - mitk::DiffusionImage::Pointer GetOutputDiffusionImage() - { - return m_OutputDiffusionImage; - } - ImageType::Pointer GetVectorImage() + ImageType::Pointer GetCorrectedDiffusionVolumes() { - return m_VectorImage; + return m_CorrectedDiffusionVolumes; } itk::Image::Pointer GetMask() { return m_MaskImage; } - //itkGetMacro(OutputDiffusionImage, mitk::DiffusionImage) - - //itkGetMacro( GradientDirectionContainer, GradientDirectionContainerType::Pointer); - protected: TensorReconstructionWithEigenvalueCorrectionFilter(); - ~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,itk::VectorImage::Pointer corrected_diffusion_temp); void CalculateAttenuation(vnl_vector org_data, vnl_vector &atten,int nof,int numberb0); - - void CalculateTensor(vnl_vector atten, vnl_vector &tensor,int nof,int numberb0); - - void CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion_temp,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min); + void CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min); void GenerateTensorImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,double what_mask,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, 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; - mitk::DiffusionImage::Pointer m_OutputDiffusionImage; - - ImageType::Pointer m_VectorImage; + ImageType::Pointer m_CorrectedDiffusionVolumes; float m_B0Threshold; - itk::Image::Pointer m_MaskImage; vnl_matrix m_PseudoInverse; vnl_matrix m_H; vnl_vector m_BVec; + + /** m_B0Mask indicates whether a volume identified by an index is B0-weighted or not */ vnl_vector m_B0Mask; - vnl_vector m_Voxdim; - int m_Flagstatus; 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 a8973cdbe7..4a5ca091c2 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.txx @@ -1,1028 +1,983 @@ /*=================================================================== 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. ===================================================================*/ -/*========================================================================= - -Program: Tensor ToolKit - TTK -Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $ -Language: C++ -Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ -Version: $Revision: 68 $ - -Copyright (c) INRIA 2010. All rights reserved. -See LICENSE.txt for details. - -This software is distributed WITHOUT ANY WARRANTY; without even -the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR -PURPOSE. See the above copyright notices for more information. - -=========================================================================*/ #ifndef _itk_TensorReconstructionWithEigenvalueCorrectionFilter_txx_ #define _itk_TensorReconstructionWithEigenvalueCorrectioFiltern_txx_ #endif #include "itkImageRegionConstIterator.h" #include -#include #include "itkImageFileWriter.h" #include "itkImage.h" #include "itkImageRegionIterator.h" namespace itk { template TensorReconstructionWithEigenvalueCorrectionFilter ::TensorReconstructionWithEigenvalueCorrectionFilter() { m_B0Threshold = 50.0; } template void TensorReconstructionWithEigenvalueCorrectionFilter ::GenerateData () { m_GradientImagePointer = static_cast< GradientImagesType * >( this->ProcessObject::GetInput(0) ); typename GradientImagesType::SizeType size = m_GradientImagePointer->GetLargestPossibleRegion().GetSize(); - int nof = m_GradientDirectionContainer->Size(); - int numberb0=0; - int cnt=0; + // 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) { numberb0++; } } - itk::Vector spacing_term = m_GradientImagePointer->GetSpacing(); - itk::Matrix direction_term = m_GradientImagePointer->GetDirection(); - - vnl_vector spacing_vnl(3); - vnl_matrix dir_vnl (3,3); - - for (int i=0;i<3;i++) - { - spacing_vnl[i]=spacing_term[i]; - - for(int j=0;j<3;j++) - { - dir_vnl[i][j]=direction_term[i][j]; - } - } - - vnl_matrix vox_dim_step (3,3); - - for (int i=0;i<3;i++) - { - for(int j=0;j<3;j++) - { - vox_dim_step[i][j]=spacing_vnl[i]*dir_vnl[i][j]; - } - } - - vnl_symmetric_eigensystem eigen_spacing(vox_dim_step); - - vnl_vector vox_dim (3); - - vox_dim[0]=eigen_spacing.get_eigenvalue(0); - vox_dim[1]=eigen_spacing.get_eigenvalue(1); - vox_dim[2]=eigen_spacing.get_eigenvalue(2); - - vox_dim=vox_dim/(vox_dim.min_value()); - - + // 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); if(vec[0]<0.0001 && vec[1]<0.0001 && vec[2]<0.0001 && vec[0]>-0.001&& vec[1]>-0.001 && vec[2]>-0.001) { + // 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); - int etbt[6] = { 0, 4, 8, 1, 5, 2 }; + //H is matrix that containes 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(); - itk::Index<3> ix; - - ImageType::Pointer corrected_diffusion = ImageType::New(); - corrected_diffusion->SetRegions(size); - corrected_diffusion->SetSpacing(m_GradientImagePointer->GetSpacing()); - corrected_diffusion->SetOrigin(m_GradientImagePointer->GetOrigin()); - corrected_diffusion->SetVectorLength(nof); - corrected_diffusion->Allocate(); - ImageType::Pointer corrected_diffusion_temp = ImageType::New(); typedef itk::VariableLengthVector VariableVectorType; VariableVectorType variableLengthVector; variableLengthVector.SetSize(nof); typedef itk::VariableLengthVector VariableVectorType; VariableVectorType corrected_single; corrected_single.SetSize(nof-1); 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->Allocate(); - - double mean_b=0.0; - + // 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(int x=0;x ix = {{x,y,z}}; + GradientVectorType pixel = m_GradientImagePointer->GetPixel(ix); for (int i=0;i m_B0Threshold) { mask->SetPixel(ix, 1); mask_cnt++; } else { mask->SetPixel(ix, 0); } } } } + + //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 ix = {{x,y,z}}; - GradientVectorType pixel2 = m_GradientImagePointer->GetPixel(ix); + GradientVectorType p = m_GradientImagePointer->GetPixel(ix); - corrected_diffusion->SetPixel(ix,pixel2); + m_CorrectedDiffusionVolumes->SetPixel(ix,p); } } } - //Sometimes the gradient voxels may contain negative values ( even if B0 voxel is > = 50 ). This must be corrected + //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. double mask_val=0.0; - vnl_vector org_vec(nof-numberb0); + vnl_vector org_vec(nof-numberb0); - int counter_corrected =0; + int counter_corrected =0; - for ( int x=0;x ix = {{x,y,z}}; - mask_val = mask->GetPixel(ix); - GradientVectorType pixel2 = corrected_diffusion->GetPixel(ix); + mask_val = mask->GetPixel(ix); + GradientVectorType pixel2 = m_CorrectedDiffusionVolumes->GetPixel(ix); - for (int i=0;i0)// we are dooing this only if the voxels are in the mask - { + if(mask_val >0) + { - for( int f=0;fSetPixel(ix, pixel2); + for (int i=0;iSetPixel(ix, pixel2); - } - } } } + } + } - typedef itk::Image< itk::DiffusionTensor3D, 3 > TensorImageType; - TensorImageType::Pointer tensorImg = TensorImageType::New(); + typename TensorImageType::Pointer tensorImg = TensorImageType::New(); tensorImg->SetRegions(m_GradientImagePointer->GetLargestPossibleRegion().GetSize()); tensorImg->SetSpacing(m_GradientImagePointer->GetSpacing()); tensorImg->SetOrigin(m_GradientImagePointer->GetOrigin()); tensorImg->Allocate(); - TensorImageType::Pointer temp_tensorImg = TensorImageType::New(); + typename TensorImageType::Pointer temp_tensorImg = TensorImageType::New(); + // Deep copy a temporary tensor image for the pre-processing methods. DeepCopyTensorImage(tensorImg,temp_tensorImg); + //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;iSetNthOutput(0, tensorImg); - m_VectorImage = corrected_diffusion; - - m_Voxdim = vox_dim; - - } template void TensorReconstructionWithEigenvalueCorrectionFilter ::SetGradientImage( GradientDirectionContainerType *gradientDirection, const GradientImagesType *gradientImage ) { - // Make sure crazy users did not call both AddGradientImage and - // SetGradientImage + 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, itk::VectorImage::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]; + int y_max=size[1]; + int z_max=size[2]; + double back_x=std::max(0,x-1); double back_y=std::max(0,y-1); double back_z=std::max(0,z-1); - int x_max=size[0];int y_max=size[1];int z_max=size[2];// converting short to int - 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);// setting constraints for neigborhood + double forth_z=std::min((z+1),z_max); double tempsum=0; itk::Index<3> ix; double temp_number=0; - double temp_mask=0;// declaration of variables - double one =1.0; - + double temp_mask=0; - for(int i=back_x; iGetPixel(ix); - //if(temp_mask > 0 && temp_mask < 2 ) - //{ - //GradientVectorType p = m_GradientImagePointer->GetPixel(ix); - GradientVectorType p = corrected_diffusion_temp->GetPixel(ix); - - double test= p[f]; - - if (test > 0.0 )// hmm this must be here becaus the method is used in multiple ocasions. Sometiems we may deal with negative values - { - if(i==x && j==y && c== z) - { - - - - } - else - { - tempsum=tempsum+p[f];// sum for calculation of mean - temp_number++;// number for calculation of mean - } + GradientVectorType p = corrected_diffusion_temp->GetPixel(ix); + double test= p[f]; + if (test > 0.0 )// taking only positive values and counting them + { + if(!(i==x && j==y && k== z)) + { + tempsum=tempsum+p[f]; + temp_number++; } - //}//end of mask condition - + } } } - }// end of size - + } + //getting back to the original position of the voxel ix[0] = x;ix[1] = y;ix[2] = z; if (temp_number <= 0.0) { - //ix[0] = x;ix[1] = y;ix[2] = z; - //GradientVectorType p = m_GradientImagePointer->GetPixel(ix); - //tempsum=p[f]; tempsum=0; mask->SetPixel(ix,0); } else { tempsum=tempsum/temp_number; } - - double ret = tempsum; - - - return ret; + 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) { 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 - void - TensorReconstructionWithEigenvalueCorrectionFilter - ::CalculateTensor(vnl_vector atten,vnl_vector &tensor, int nof,int numberb0) - { - for (int i=0;i double TensorReconstructionWithEigenvalueCorrectionFilter ::CheckNegatives ( itk::Size<3> size, itk::Image::Pointer mask, 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; double pixel=0; itk::Index<3> ix; itk::DiffusionTensor3D ten; vnl_matrix temp_tensor(3,3); - vnl_vector eigen_vals(3);// declaration of important structures and variables + vnl_vector eigen_vals(3); vnl_vector tensor (6); - + // for every pixel from the image for (int x=0;xGetPixel(ix); + // but only if previously marked as bad one-negative eigen value + if(pixel > 1) - {// but only if previously marked as bad or potentially bad one + { ten = tensorImg->GetPixel(ix); + // changing order from tensor structure in MITK into vnl structure 3x3 symmetric tensor matrix + 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); 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]; + //checking negativity of tensor eigenvalues - vnl_symmetric_eigensystem eigen_tensor(temp_tensor);//creating an eigensystem out of 3x3 symmetric tensor matrix + 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);// calculating eigenvalues for current tensor + eigen_vals[2]=eigen_tensor.get_eigenvalue(2); + + //comparison to 0.01 instead of 0 was proposed by O.Pasternak if( eigen_vals[0]>0.01 && eigen_vals[1]>0.01 && eigen_vals[2]>0.01) { mask->SetPixel(ix,1); } else { - badvoxels++;// increasing a number of detected bad + badvoxels++; } } } } } - double ret = badvoxels;//returning just the number of bad voxels + double ret = badvoxels; return ret; - }// end of void check negativity + } template void TensorReconstructionWithEigenvalueCorrectionFilter - ::CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion_temp,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min) + ::CorrectDiffusionImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,vnl_vector< double> pixel_max,vnl_vector< double> pixel_min) { - // in this method the diffusion image temporray is corrected with use of information from updated mask. The diffusion image is tempora - //-ry while we can't use corrected voxel to correct other voxel in the same iteration. + // 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 itk::Index<3> ix; vnl_vector org_data(nof-numberb0); vnl_vector atten(nof-numberb0); - double cnt_atten=0;// declaration of important variables + double cnt_atten=0; for (int z=0;zGetPixel(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]) { - org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion_temp); + org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion); } cnt_atten++; } + //smoothing B0 if(m_B0Mask[f]==1) { - - org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion_temp); - + org_data[f] = CheckNeighbours(x,y,z,f,size,mask,corrected_diffusion); } - }//end for - - + } for (int i=0;iSetPixel(ix, pt); + corrected_diffusion->SetPixel(ix, pt); } else { GradientVectorType pt = corrected_diffusion->GetPixel(ix); - corrected_diffusion_temp->SetPixel(ix, pt); + corrected_diffusion->SetPixel(ix, pt); } } } } - std::cout << "break point"; + } template void TensorReconstructionWithEigenvalueCorrectionFilter ::GenerateTensorImage(int nof,int numberb0,itk::Size<3> size,itk::VectorImage::Pointer corrected_diffusion,itk::Image::Pointer mask,double what_mask,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-numberb0); vnl_vector atten(nof-numberb0); vnl_vector tensor(6); itk::DiffusionTensor3D ten; double mask_val=0; for (int x=0;xGetPixel(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) { 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;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; 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) { - // in this method voxels in the mask that poses a certain value are substituded by other value. It is connected to the - //idea in the algorithm. It is called more tna once so separate method spares some lines of code. + // The method changes voxels in the mask that poses a certain value with other value. + itk::Index<3> ix; + double temp_mask_value=0; - itk::Index<3> ix; - double temp_mask_value=0; - - for(int x=0;xGetPixel(ix); if(temp_mask_value>previous_mask) { - mask->SetPixel(ix,set_mask); + mask->SetPixel(ix,set_mask); } + } } } - } -}//end of turn mas + } 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()); inputIterator.GoToBegin(); outputIterator.GoToBegin(); 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) { temp_tensorImg->SetSpacing(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()) { outputIterator.Set(inputIterator.Get()); ++inputIterator; ++outputIterator; } } } // end of namespace