diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp index 4410ee0611..2a22ba2f3a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp @@ -1,318 +1,324 @@ /*=================================================================== 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 MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP #define MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP #include "mitkDWIHeadMotionCorrectionFilter.h" #include "itkSplitDWImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkPyramidImageRegistrationMethod.h" +//#include "mitkRegistrationMethodITK4.h" #include "mitkImageToDiffusionImageSource.h" #include "mitkDiffusionImageCorrectionFilter.h" #include #include #include "mitkIOUtil.h" #include template< typename DiffusionPixelType> mitk::DWIHeadMotionCorrectionFilter - ::DWIHeadMotionCorrectionFilter(): - m_CurrentStep(0), - m_Steps(100), - m_IsInValidState(true), - m_AbortRegistration(false) + ::DWIHeadMotionCorrectionFilter() + : m_CurrentStep(0), + m_Steps(100), + m_IsInValidState(true), + m_AbortRegistration(false), + m_AverageUnweighted(true) { } template< typename DiffusionPixelType> void mitk::DWIHeadMotionCorrectionFilter - ::GenerateData() +::GenerateData() { typedef itk::SplitDWImageFilter< DiffusionPixelType, DiffusionPixelType> SplitFilterType; DiffusionImageType* input = const_cast(this->GetInput(0)); + typedef mitk::PyramidImageRegistrationMethod RegistrationMethod; + // typedef mitk::RegistrationMethodITKV4 RegistrationMethod + unsigned int numberOfSteps = input->GetVectorImage()->GetNumberOfComponentsPerPixel () ; m_Steps = numberOfSteps; // // (1) Extract the b-zero images to a 3d+t image, register them by NCorr metric and // rigid registration : they will then be used are reference image for registering // the gradient images // typedef itk::B0ImageExtractionToSeparateImageFilter< DiffusionPixelType, DiffusionPixelType> B0ExtractorType; typename B0ExtractorType::Pointer b0_extractor = B0ExtractorType::New(); b0_extractor->SetInput( input->GetVectorImage() ); b0_extractor->SetDirections( input->GetDirections() ); b0_extractor->Update(); mitk::Image::Pointer b0Image = mitk::Image::New(); b0Image->InitializeByItk( b0_extractor->GetOutput() ); b0Image->SetImportChannel( b0_extractor->GetOutput()->GetBufferPointer(), - mitk::Image::CopyMemory ); + mitk::Image::CopyMemory ); // (2.1) Use the extractor to access the extracted b0 volumes mitk::ImageTimeSelector::Pointer t_selector = - mitk::ImageTimeSelector::New(); + mitk::ImageTimeSelector::New(); t_selector->SetInput( b0Image ); t_selector->SetTimeNr(0); t_selector->Update(); // first unweighted image as reference space for the registration mitk::Image::Pointer b0referenceImage = t_selector->GetOutput(); + const unsigned int numberOfb0Images = b0Image->GetTimeSteps(); - mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New(); - registrationMethod->SetFixedImage( b0referenceImage ); - registrationMethod->SetTransformToRigid(); + // register b-zeros only if the flag to average is set, use the first one otherwise + if( m_AverageUnweighted && numberOfb0Images > 1) + { - // the unweighted images are of same modality - registrationMethod->SetCrossModalityOff(); + RegistrationMethod::Pointer registrationMethod = RegistrationMethod::New(); + registrationMethod->SetFixedImage( b0referenceImage ); + registrationMethod->SetTransformToRigid(); - // use the advanced (windowed sinc) interpolation - registrationMethod->SetUseAdvancedInterpolation(true); + // the unweighted images are of same modality + registrationMethod->SetCrossModalityOff(); - // Initialize the temporary output image - mitk::Image::Pointer registeredB0Image = b0Image->Clone(); - const unsigned int numberOfb0Images = b0Image->GetTimeSteps(); + // use the advanced (windowed sinc) interpolation + registrationMethod->SetUseAdvancedInterpolation(true); + + // Initialize the temporary output image + mitk::Image::Pointer registeredB0Image = b0Image->Clone(); - if( numberOfb0Images > 1) - { mitk::ImageTimeSelector::Pointer t_selector2 = - mitk::ImageTimeSelector::New(); + mitk::ImageTimeSelector::New(); t_selector2->SetInput( b0Image ); for( unsigned int i=1; iSetTimeNr(i); - t_selector2->Update(); + t_selector2->SetTimeNr(i); + t_selector2->Update(); - registrationMethod->SetMovingImage( t_selector2->GetOutput() ); + registrationMethod->SetMovingImage( t_selector2->GetOutput() ); - try - { - MITK_INFO << " === (" << i <<"/"<< numberOfb0Images-1 << ") :: Starting registration"; - registrationMethod->Update(); - } - catch( const itk::ExceptionObject& e) - { - m_IsInValidState = false; - mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); - } + try + { + MITK_INFO << " === (" << i <<"/"<< numberOfb0Images-1 << ") :: Starting registration"; + registrationMethod->Update(); + } + catch( const itk::ExceptionObject& e) + { + m_IsInValidState = false; + mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); + } // import volume to the inter-results mitk::ImageWriteAccessor imac(registrationMethod->GetResampledMovingImage()); registeredB0Image->SetImportVolume( imac.GetData(), i, 0, mitk::Image::ReferenceMemory ); - - } + } // use the accumulateImageFilter as provided by the ItkAccumulateFilter method in the header file AccessFixedDimensionByItk_1(registeredB0Image, ItkAccumulateFilter, (4), b0referenceImage ); } // // (2) Split the diffusion image into a 3d+t regular image, extract only the weighted images // typename SplitFilterType::Pointer split_filter = SplitFilterType::New(); split_filter->SetInput (input->GetVectorImage() ); - split_filter->SetExtractAllAboveThreshold(20, input->GetBValueMap() ); + split_filter->SetExtractAllAboveThreshold(8, input->GetB_ValueMap() ); try { split_filter->Update(); } catch( const itk::ExceptionObject &e) { m_IsInValidState = false; mitkThrow() << " Caught exception from SplitImageFilter : " << e.what(); } mitk::Image::Pointer splittedImage = mitk::Image::New(); splittedImage->InitializeByItk( split_filter->GetOutput() ); splittedImage->SetImportChannel( split_filter->GetOutput()->GetBufferPointer(), - mitk::Image::CopyMemory ); + mitk::Image::CopyMemory ); // // (3) Use again the time-selector to access the components separately in order // to perform the registration of Image -> unweighted reference // - mitk::PyramidImageRegistrationMethod::Pointer weightedRegistrationMethod - = mitk::PyramidImageRegistrationMethod::New(); + RegistrationMethod::Pointer weightedRegistrationMethod + = RegistrationMethod::New(); weightedRegistrationMethod->SetTransformToAffine(); weightedRegistrationMethod->SetCrossModalityOn(); // // - (3.1) Set the reference image // - a single b0 image // - average over the registered b0 images if multiple present // weightedRegistrationMethod->SetFixedImage( b0referenceImage ); // use the advanced (windowed sinc) interpolation weightedRegistrationMethod->SetUseAdvancedInterpolation(true); // // - (3.2) Register all timesteps in the splitted image onto the first reference // unsigned int maxImageIdx = splittedImage->GetTimeSteps(); mitk::TimeGeometry* tsg = splittedImage->GetTimeGeometry(); mitk::ProportionalTimeGeometry* ptg = dynamic_cast(tsg); ptg->Expand(maxImageIdx+1); ptg->SetTimeStepGeometry( ptg->GetGeometryForTimeStep(0), maxImageIdx ); mitk::Image::Pointer registeredWeighted = mitk::Image::New(); registeredWeighted->Initialize( splittedImage->GetPixelType(0), *tsg ); // insert the first unweighted reference as the first volume mitk::ImageWriteAccessor imac(b0referenceImage); registeredWeighted->SetImportVolume( imac.GetData(), 0,0, mitk::Image::CopyMemory ); // mitk::Image::Pointer registeredWeighted = splittedImage->Clone(); // this time start at 0, we have only gradient images in the 3d+t file // the reference image comes form an other image mitk::ImageTimeSelector::Pointer t_selector_w = - mitk::ImageTimeSelector::New(); + mitk::ImageTimeSelector::New(); t_selector_w->SetInput( splittedImage ); // store the rotation parts of the transformations in a vector - typedef mitk::PyramidImageRegistrationMethod::TransformMatrixType MatrixType; + typedef RegistrationMethod::TransformMatrixType MatrixType; std::vector< MatrixType > estimated_transforms; for( unsigned int i=0; iSetTimeNr(i); t_selector_w->Update(); weightedRegistrationMethod->SetMovingImage( t_selector_w->GetOutput() ); try { MITK_INFO << " === (" << i+1 <<"/"<< maxImageIdx << ") :: Starting registration"; weightedRegistrationMethod->Update(); } catch( const itk::ExceptionObject& e) { m_IsInValidState = false; mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // allow expansion mitk::ImageWriteAccessor imac(weightedRegistrationMethod->GetResampledMovingImage()); registeredWeighted->SetImportVolume( imac.GetData(), i+1, 0, mitk::Image::CopyMemory); estimated_transforms.push_back( weightedRegistrationMethod->GetLastRotationMatrix() ); } // // (4) Cast the resulting image back to an diffusion weighted image // typename DiffusionImageType::GradientDirectionContainerType *gradients = input->GetDirections(); typename DiffusionImageType::GradientDirectionContainerType::Pointer gradients_new = - DiffusionImageType::GradientDirectionContainerType::New(); + DiffusionImageType::GradientDirectionContainerType::New(); typename DiffusionImageType::GradientDirectionType bzero_vector; bzero_vector.fill(0); // compose the direction vector // - no direction for the first image // - correct ordering of the directions based on the index list gradients_new->push_back( bzero_vector ); typename SplitFilterType::IndexListType index_list = split_filter->GetIndexList(); typename SplitFilterType::IndexListType::const_iterator lIter = index_list.begin(); while( lIter != index_list.end() ) { gradients_new->push_back( gradients->at( *lIter ) ); ++lIter; } typename mitk::ImageToDiffusionImageSource< DiffusionPixelType >::Pointer caster = - mitk::ImageToDiffusionImageSource< DiffusionPixelType >::New(); + mitk::ImageToDiffusionImageSource< DiffusionPixelType >::New(); caster->SetImage( registeredWeighted ); caster->SetBValue( input->GetReferenceBValue() ); caster->SetGradientDirections( gradients_new.GetPointer() ); try { caster->Update(); } catch( const itk::ExceptionObject& e) { m_IsInValidState = false; MITK_ERROR << "Casting back to diffusion image failed: "; mitkThrow() << "Subprocess failed with exception: " << e.what(); } // // (5) Adapt the gradient directions according to the estimated transforms // typedef mitk::DiffusionImageCorrectionFilter< DiffusionPixelType > CorrectionFilterType; typename CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); OutputImagePointerType output = caster->GetOutput(); corrector->SetImage( output ); corrector->CorrectDirections( estimated_transforms ); // // (6) Pass the corrected image to the filters output port // m_CurrentStep += 1; this->GetOutput()->SetVectorImage(output->GetVectorImage()); this->GetOutput()->SetReferenceBValue(output->GetReferenceBValue()); this->GetOutput()->SetMeasurementFrame(output->GetMeasurementFrame()); this->GetOutput()->SetDirections(output->GetDirections()); this->GetOutput()->InitializeFromVectorImage(); this->GetOutput()->Modified(); } #endif // MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h index 888c5719f4..d8065ff27d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h @@ -1,130 +1,132 @@ /*=================================================================== 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 MITKDWIHEADMOTIONCORRECTIONFILTER_H #define MITKDWIHEADMOTIONCORRECTIONFILTER_H #include "mitkDiffusionImageToDiffusionImageFilter.h" #include #include #include "mitkITKImageImport.h" namespace mitk { /** * @class DWIHeadMotionCorrectionFilter * * @brief Performs standard head-motion correction by using affine registration of the gradient images. * * (Head) motion correction is a essential pre-processing step before performing any further analysis of a diffusion-weighted * images since all model fits ( tensor, QBI ) rely on an aligned diffusion-weighted dataset. The correction is done in two steps. First the * unweighted images ( if multiple present ) are separately registered on the first one by means of rigid registration and normalized correlation * as error metric. Second, the weighted gradient images are registered to the unweighted reference ( computed as average from the aligned images from first step ) * by an affine transformation using the MattesMutualInformation metric as optimizer guidance. * */ template< typename DiffusionPixelType> class DWIHeadMotionCorrectionFilter : public DiffusionImageToDiffusionImageFilter< DiffusionPixelType > { public: // class macros mitkClassMacro( DWIHeadMotionCorrectionFilter, DiffusionImageToDiffusionImageFilter ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkGetMacro( Steps, unsigned long ) itkGetMacro( CurrentStep, unsigned long ) itkGetMacro( IsInValidState, bool) itkSetMacro( AbortRegistration, bool ) + itkSetMacro( AverageUnweighted, bool ) // public typedefs typedef typename Superclass::InputImageType DiffusionImageType; typedef typename Superclass::InputImagePointerType DiffusionImagePointerType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImagePointerType OutputImagePointerType; protected: DWIHeadMotionCorrectionFilter(); virtual ~DWIHeadMotionCorrectionFilter() {} virtual void GenerateData(); unsigned long m_CurrentStep; unsigned long m_Steps; bool m_IsInValidState; ///< Whether the filter is in a valid state, false if error occured bool m_AbortRegistration; ///< set flag to abort + bool m_AverageUnweighted; - /** - * @brief Averages an 3d+t image along the time axis. - * - * The method uses the AccumulateImageFilter as provided by ITK and collapses the given 3d+t image - * to an 3d image while computing the average over the time axis for each of the spatial voxels. - */ - template< typename TPixel, unsigned int VDimensions> - static void ItkAccumulateFilter( - const itk::Image< TPixel, VDimensions>* image, - mitk::Image::Pointer& output) +}; + +/** + * @brief Averages an 3d+t image along the time axis. + * + * The method uses the AccumulateImageFilter as provided by ITK and collapses the given 3d+t image + * to an 3d image while computing the average over the time axis for each of the spatial voxels. + */ +template< typename TPixel, unsigned int VDimensions> +static void ItkAccumulateFilter( + const itk::Image< TPixel, VDimensions>* image, + mitk::Image::Pointer& output) +{ + // input 3d+t --> output 3d + typedef itk::Image< TPixel, 4> InputItkType; + typedef itk::Image< TPixel, 3> OutputItkType; + // the accumulate filter requires the same dimension in output and input image + typedef typename itk::AccumulateImageFilter< InputItkType, InputItkType > FilterType; + + typename FilterType::Pointer filter = FilterType::New(); + filter->SetInput( image ); + filter->SetAccumulateDimension( 3 ); + filter->SetAverage( true ); + + // we need to extract the volume to reduce the size from 4 to 3 for further processing + typedef typename itk::ExtractImageFilter< InputItkType, OutputItkType > ExtractFilterType; + typename ExtractFilterType::Pointer extractor = ExtractFilterType::New(); + extractor->SetInput( filter->GetOutput() ); + extractor->SetDirectionCollapseToIdentity(); + + typename InputItkType::RegionType extractRegion = image->GetLargestPossibleRegion(); + // crop out the time axis + extractRegion.SetSize( 3, 0); + extractor->SetExtractionRegion( extractRegion ); + + try + { + extractor->Update(); + } + catch( const itk::ExceptionObject& e) { - // input 3d+t --> output 3d - typedef itk::Image< TPixel, 4> InputItkType; - typedef itk::Image< TPixel, 3> OutputItkType; - // the accumulate filter requires the same dimension in output and input image - typedef typename itk::AccumulateImageFilter< InputItkType, InputItkType > FilterType; - - typename FilterType::Pointer filter = FilterType::New(); - filter->SetInput( image ); - filter->SetAccumulateDimension( 3 ); - filter->SetAverage( true ); - - // we need to extract the volume to reduce the size from 4 to 3 for further processing - typedef typename itk::ExtractImageFilter< InputItkType, OutputItkType > ExtractFilterType; - typename ExtractFilterType::Pointer extractor = ExtractFilterType::New(); - extractor->SetInput( filter->GetOutput() ); - extractor->SetDirectionCollapseToIdentity(); - - typename InputItkType::RegionType extractRegion = image->GetLargestPossibleRegion(); - // crop out the time axis - extractRegion.SetSize( 3, 0); - extractor->SetExtractionRegion( extractRegion ); - - try - { - extractor->Update(); - } - catch( const itk::ExceptionObject& e) - { - mitkThrow() << " Exception while averaging: " << e.what(); - } - - output = mitk::GrabItkImageMemory( extractor->GetOutput() ); + mitkThrow() << " Exception while averaging: " << e.what(); } -}; + output = mitk::GrabItkImageMemory( extractor->GetOutput() ); +} } //end namespace mitk #include "mitkDWIHeadMotionCorrectionFilter.cpp" #endif // MITKDWIHEADMOTIONCORRECTIONFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp index 418a5a1bd7..ad71e6b9a5 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp @@ -1,180 +1,181 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPyramidImageRegistrationMethod.h" #include "mitkException.h" #include "mitkImageAccessByItk.h" mitk::PyramidImageRegistrationMethod::PyramidImageRegistrationMethod() : m_FixedImage(NULL), m_MovingImage(NULL), m_CrossModalityRegistration(true), m_UseAffineTransform(true), m_UseWindowedSincInterpolator(false), m_UseNearestNeighborInterpolator(false), m_UseMask(false), m_EstimatedParameters(NULL), - m_Verbose(false), - m_InitializeByGeometry(false) + m_InitialParameters(NULL), + m_Verbose(false) { } mitk::PyramidImageRegistrationMethod::~PyramidImageRegistrationMethod() { if( m_EstimatedParameters != NULL) { delete [] m_EstimatedParameters; } + } void mitk::PyramidImageRegistrationMethod::SetFixedImage(mitk::Image::Pointer fixed) { if( fixed.IsNotNull() ) { m_FixedImage = fixed; } } void mitk::PyramidImageRegistrationMethod::SetMovingImage(mitk::Image::Pointer moving) { if( moving.IsNotNull() ) { m_MovingImage = moving; } } void mitk::PyramidImageRegistrationMethod::SetFixedImageMask(mitk::Image::Pointer mask) { m_FixedImageMask = mask; } void mitk::PyramidImageRegistrationMethod::Update() { if( m_MovingImage.IsNull() ) { mitkThrow() << " Moving image is null"; } if( m_FixedImage.IsNull() ) { mitkThrow() << " Moving image is null"; } unsigned int allowedDimension = 3; if( m_FixedImage->GetDimension() != allowedDimension || m_MovingImage->GetDimension() != allowedDimension ) { mitkThrow() << " Only 3D Images supported."; } // // One possibility: use the FixedDimesnionByItk, but this instantiates ALL possible // pixel type combinations! // AccessTwoImagesFixedDimensionByItk( m_FixedImage, m_MovingImage, RegisterTwoImages, 3); // in helper: TypeSubset : short, float AccessTwoImagesFixedDimensionTypeSubsetByItk( m_FixedImage, m_MovingImage, RegisterTwoImages, 3); } mitk::PyramidImageRegistrationMethod::TransformMatrixType mitk::PyramidImageRegistrationMethod ::GetLastRotationMatrix() { TransformMatrixType output; if( m_EstimatedParameters == NULL ) { output.set_identity(); return output; } typedef itk::MatrixOffsetTransformBase< double, 3, 3> BaseTransformType; BaseTransformType::Pointer base_transform = BaseTransformType::New(); if( this->m_UseAffineTransform ) { typedef itk::AffineTransform< double > TransformType; TransformType::Pointer transform = TransformType::New(); TransformType::ParametersType affine_params( TransformType::ParametersDimension ); this->GetParameters( &affine_params[0] ); transform->SetParameters( affine_params ); base_transform = transform; } else { typedef itk::Euler3DTransform< double > RigidTransformType; RigidTransformType::Pointer rtransform = RigidTransformType::New(); RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); this->GetParameters( &rigid_params[0] ); rtransform->SetParameters( rigid_params ); base_transform = rtransform; } return base_transform->GetMatrix().GetVnlMatrix(); } mitk::Image::Pointer mitk::PyramidImageRegistrationMethod ::GetResampledMovingImage() { mitk::Image::Pointer output = mitk::Image::New(); //output->Initialize( this->m_FixedImage ); AccessFixedDimensionByItk_1( this->m_MovingImage, ResampleMitkImage, 3, output ); return output; } mitk::Image::Pointer mitk::PyramidImageRegistrationMethod::GetResampledMovingImage(mitk::Image::Pointer movingImage, double* transform) { mitk::Image::Pointer output = mitk::Image::New(); unsigned int dim = 12; if( !m_UseAffineTransform ) dim = 6; if (m_EstimatedParameters == NULL) m_EstimatedParameters = new double[dim]; double tmpParams[12]; // save and set temporal transformation values for( unsigned int i=0; i #include #include "itkImageMaskSpatialObject.h" #include "itkNotImageFilter.h" #include #include #include "mitkImage.h" #include "mitkBaseProcess.h" #include "mitkPyramidRegistrationMethodHelper.h" #include #include "mitkImageToItk.h" #include "mitkImageCast.h" #include "mitkImageCaster.h" #include "mitkITKImageImport.h" #include "mitkIOUtil.h" namespace mitk { /** * @brief The PyramidImageRegistration class implements a multi-scale registration method * * The PyramidImageRegistration class is suitable for aligning (f.e.) brain MR images. The method offers two * transform types * - Rigid: optimizing translation and rotation only and * - Affine ( default ): with scaling in addition ( 12 DOF ) * * The error metric is internally chosen based on the selected task type ( @sa SetCrossModalityOn ) * It uses * - MattesMutualInformation for CrossModality=on ( default ) and * - NormalizedCorrelation for CrossModality=off. */ class MitkDiffusionCore_EXPORT PyramidImageRegistrationMethod : public itk::Object { public: /** Typedefs */ mitkClassMacro(PyramidImageRegistrationMethod, itk::Object) /** Smart pointer support */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) + typedef itk::OptimizerParameters ParametersType; + /** Typedef for the transformation matrix, corresponds to the InternalMatrixType from ITK transforms */ typedef vnl_matrix_fixed< double, 3, 3> TransformMatrixType; typedef itk::AffineTransform< double > AffineTransformType; typedef itk::Euler3DTransform< double > RigidTransformType; /** Registration is between modalities - will take MattesMutualInformation as error metric */ void SetCrossModalityOn() { m_CrossModalityRegistration = true; } /** Registration is between modalities - will take NormalizedCorrelation as error metric */ void SetCrossModalityOff() { m_CrossModalityRegistration = false; } /** Turn the cross-modality on/off */ void SetCrossModality(bool flag) { if( flag ) this->SetCrossModalityOn(); else this->SetCrossModalityOff(); } /** Controll the verbosity of the registration output * * for true, each iteration step of the optimizer is watched and passed to the std::cout */ void SetVerbose(bool flag) { if( flag ) this->SetVerboseOn(); else this->SetVerboseOff(); } /** Turn verbosity on, \sa SetVerbose */ void SetVerboseOn() { m_Verbose = true; } /** Turn verbosity off, \sa SetVerbose */ void SetVerboseOff() { m_Verbose = false; } /** A rigid ( 6dof : translation + rotation ) transform is optimized */ void SetTransformToRigid() { m_UseAffineTransform = false; } /** An affine ( 12dof : Rigid + scale + skew ) transform is optimized */ void SetTransformToAffine() { m_UseAffineTransform = true; } /** Input image, the reference one */ void SetFixedImage( mitk::Image::Pointer ); /** Input image, the one to be transformed */ void SetMovingImage( mitk::Image::Pointer ); /** Fixed image mask, excludes the masked voxels from the registration metric*/ void SetFixedImageMask( mitk::Image::Pointer mask); void SetInitializeByGeometry(bool flag) { m_InitializeByGeometry = flag; } void Update(); /** * @brief Get the number of parameters optimized ( 12 or 6 ) * * @return number of paramters */ unsigned int GetNumberOfParameters() { unsigned int retValue = 12; if(!m_UseAffineTransform) retValue = 6; return retValue; } /** * @brief Copies the estimated parameters to the given array * @param paramArray, target array for copy, make sure to allocate enough space */ void GetParameters( double* paramArray) { if( m_EstimatedParameters == NULL ) { mitkThrow() << "No parameters were estimated yet, call Update() first."; } unsigned int dim = 12; if( !m_UseAffineTransform ) dim = 6; for( unsigned int i=0; i * * It returns identity if the internal parameters are not-yet allocated ( i.e. the filter did not run yet ) * * @return \sa TransformMatrixType */ TransformMatrixType GetLastRotationMatrix(); protected: PyramidImageRegistrationMethod(); ~PyramidImageRegistrationMethod(); /** Fixed image, used as reference for registration */ mitk::Image::Pointer m_FixedImage; /** Moving image, will be transformed */ mitk::Image::Pointer m_MovingImage; mitk::Image::Pointer m_FixedImageMask; bool m_CrossModalityRegistration; bool m_UseAffineTransform; bool m_UseWindowedSincInterpolator; bool m_UseNearestNeighborInterpolator; bool m_UseMask; double* m_EstimatedParameters; + ParametersType m_InitialParameters; + /** Control the verbosity of the regsitistration output */ bool m_Verbose; bool m_InitializeByGeometry; template void RegisterTwoImages(itk::Image* itkImage1, itk::Image* itkImage2) { typedef typename itk::Image< TPixel1, VImageDimension1> FixedImageType; typedef typename itk::Image< TPixel2, VImageDimension2> MovingImageType; typedef typename itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType; typedef itk::RegularStepGradientDescentOptimizer OptimizerType; typedef itk::MatrixOffsetTransformBase< double, VImageDimension1, VImageDimension2 > BaseTransformType; typedef typename itk::MattesMutualInformationImageToImageMetric< FixedImageType, MovingImageType > MMIMetricType; typedef typename itk::NormalizedCorrelationImageToImageMetric< FixedImageType, MovingImageType > NCMetricType; typedef typename itk::ImageToImageMetric< FixedImageType, MovingImageType> BaseMetricType; typename itk::LinearInterpolateImageFunction::Pointer interpolator = itk::LinearInterpolateImageFunction::New(); typename BaseMetricType::Pointer metric; + unsigned int glob_max_threads = itk::MultiThreader::GetGlobalMaximumNumberOfThreads(); + itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1); + if( m_CrossModalityRegistration ) { metric = MMIMetricType::New(); + //metric->SetNumberOfThreads( 6 ); } else { metric = NCMetricType::New(); } // initial parameter dimension ( affine = default ) unsigned int paramDim = 12; typename BaseTransformType::Pointer transform; if( m_UseAffineTransform ) { transform = AffineTransformType::New(); } else { transform = RigidTransformType::New(); paramDim = 6; } typename BaseTransformType::ParametersType initialParams(paramDim); initialParams.Fill(0); if( m_UseAffineTransform ) { initialParams[0] = initialParams[4] = initialParams[8] = 1; } + +/* FIXME : Review removal, occured during rebasing on master <30 Jan> + if( m_InitialParameters.Size() == paramDim ) + { + initialParams = m_InitialParameters; + } +*/ typename FixedImageType::Pointer referenceImage = itkImage1; typename MovingImageType::Pointer movingImage = itkImage2; typename FixedImageType::RegionType maskedRegion = referenceImage->GetLargestPossibleRegion(); typename RegistrationType::Pointer registration = RegistrationType::New(); unsigned int max_pyramid_lvl = 3; unsigned int max_schedule_val = 4; typename FixedImageType::RegionType::SizeType image_size = referenceImage->GetLargestPossibleRegion().GetSize(); unsigned int min_value = std::min( image_size[0], std::min( image_size[1], image_size[2])); // condition for the top level pyramid image float optmaxstep = 6; float optminstep = 0.05f; unsigned int iterations = 100; if( min_value / max_schedule_val < 12 ) { max_schedule_val /= 2; + max_pyramid_lvl--; optmaxstep *= 0.25f; optminstep *= 0.1f; } typename RegistrationType::ScheduleType fixedSchedule(max_pyramid_lvl,3); fixedSchedule.Fill(1); fixedSchedule[0][0] = max_schedule_val; fixedSchedule[0][1] = max_schedule_val; fixedSchedule[0][2] = max_schedule_val; for( unsigned int i=1; iSetScales( optScales ); optimizer->SetInitialPosition( initialParams ); optimizer->SetNumberOfIterations( iterations ); optimizer->SetGradientMagnitudeTolerance( 1e-4 ); optimizer->SetRelaxationFactor( 0.7 ); optimizer->SetMaximumStepLength( optmaxstep ); optimizer->SetMinimumStepLength( optminstep ); // Initializing by Geometry if (!m_UseAffineTransform && m_InitializeByGeometry) { typedef itk::CenteredVersorTransformInitializer< FixedImageType, MovingImageType> TransformInitializerType; typename TransformInitializerType::TransformType::Pointer rigidTransform = TransformInitializerType::TransformType::New() ; MITK_INFO << "Initializer starting at : " << rigidTransform->GetParameters(); typename TransformInitializerType::Pointer initializer = TransformInitializerType::New(); initializer->SetTransform( rigidTransform); initializer->SetFixedImage( referenceImage.GetPointer() ); initializer->SetMovingImage( movingImage.GetPointer() ); initializer->MomentsOn(); initializer->InitializeTransform(); MITK_INFO << "Initialized Rigid position : " << rigidTransform->GetParameters(); initialParams[3]=rigidTransform->GetParameters()[3]; initialParams[4]=rigidTransform->GetParameters()[4]; initialParams[5]=rigidTransform->GetParameters()[5]; } // add observer tag if verbose */ unsigned long vopt_tag = 0; if(m_Verbose) { MITK_INFO << " :: Starting at " << initialParams; MITK_INFO << " :: Scales = " << optScales; OptimizerIterationCommand< OptimizerType >::Pointer iterationObserver = OptimizerIterationCommand::New(); vopt_tag = optimizer->AddObserver( itk::IterationEvent(), iterationObserver ); } // INPUT IMAGES registration->SetFixedImage( referenceImage ); registration->SetFixedImageRegion( maskedRegion ); registration->SetMovingImage( movingImage ); registration->SetSchedules( fixedSchedule, fixedSchedule); // SET MASKED AREA typedef itk::Image BinaryImageType; BinaryImageType::Pointer itkFixedImageMask = BinaryImageType::New(); itk::ImageMaskSpatialObject< 3 >::Pointer fixedMaskSpatialObject = itk::ImageMaskSpatialObject< 3 >::New(); if (m_UseMask) { CastToItkImage(m_FixedImageMask,itkFixedImageMask); itk::NotImageFilter::Pointer notFilter = itk::NotImageFilter::New(); notFilter->SetInput(itkFixedImageMask); notFilter->Update(); fixedMaskSpatialObject->SetImage( notFilter->GetOutput() ); metric->SetFixedImageMask( fixedMaskSpatialObject ); } // OTHER INPUTS registration->SetMetric( metric ); registration->SetOptimizer( optimizer ); registration->SetTransform( transform.GetPointer() ); registration->SetInterpolator( interpolator ); registration->SetInitialTransformParameters( initialParams ); typename PyramidOptControlCommand::Pointer pyramid_observer = PyramidOptControlCommand::New(); - registration->AddObserver( itk::IterationEvent(), pyramid_observer ); + + unsigned int pyramid_tag = 0; + if(m_Verbose) + { + pyramid_tag = registration->AddObserver( itk::IterationEvent(), pyramid_observer ); + } try { registration->Update(); } catch (itk::ExceptionObject &e) { + registration->Print( std::cout ); + /*registration->Print( std::cout ); + MITK_INFO << "========== reference "; + itkImage1->Print( std::cout ); + + MITK_INFO << "========== moving "; + itkImage2->Print( std::cout );*/ + MITK_ERROR << "[Registration Update] Caught ITK exception: "; mitkThrow() << "Registration failed with exception: " << e.what(); } if( m_EstimatedParameters != NULL) { delete [] m_EstimatedParameters; } m_EstimatedParameters = new double[paramDim]; typename BaseTransformType::ParametersType finalParameters = registration->GetLastTransformParameters(); for( unsigned int i=0; iGetValue() << " :: " << finalParameters; + // remove optimizer tag if used if( vopt_tag ) { optimizer->RemoveObserver( vopt_tag ); } + if( pyramid_tag ) + { + registration->RemoveObserver( pyramid_tag ); + } + + itk::MultiThreader::SetGlobalMaximumNumberOfThreads( glob_max_threads ); } template< typename TPixel, unsigned int VDimension> void ResampleMitkImage( itk::Image* itkImage, mitk::Image::Pointer& outputImage ) { typedef typename itk::Image< TPixel, VDimension> ImageType; typedef typename itk::ResampleImageFilter< ImageType, ImageType, double> ResampleImageFilterType; typedef itk::LinearInterpolateImageFunction< ImageType, double > InterpolatorType; typename InterpolatorType::Pointer linear_interpolator = InterpolatorType::New(); typedef itk::NearestNeighborInterpolateImageFunction< ImageType, double > NearestNeighborInterpolatorType; typename NearestNeighborInterpolatorType::Pointer nn_interpolator = NearestNeighborInterpolatorType::New(); typedef itk::WindowedSincInterpolateImageFunction< ImageType, 7> WindowedSincInterpolatorType; typename WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New(); typename ImageType::Pointer reference_image = ImageType::New(); CastToItkImage(m_FixedImage,reference_image); typedef itk::MatrixOffsetTransformBase< double, 3, 3> BaseTransformType; BaseTransformType::Pointer base_transform = BaseTransformType::New(); if( this->m_UseAffineTransform ) { typedef itk::AffineTransform< double > TransformType; TransformType::Pointer transform = TransformType::New(); TransformType::ParametersType affine_params( TransformType::ParametersDimension ); this->GetParameters( &affine_params[0] ); transform->SetParameters( affine_params ); base_transform = transform; } // Rigid else { typedef itk::Euler3DTransform< double > RigidTransformType; RigidTransformType::Pointer rtransform = RigidTransformType::New(); RigidTransformType::ParametersType rigid_params( RigidTransformType::ParametersDimension ); this->GetParameters( &rigid_params[0] ); rtransform->SetParameters( rigid_params ); base_transform = rtransform; } typename ResampleImageFilterType::Pointer resampler = ResampleImageFilterType::New(); resampler->SetInterpolator( linear_interpolator ); if( m_UseWindowedSincInterpolator ) resampler->SetInterpolator( sinc_interpolator ); if ( m_UseNearestNeighborInterpolator) resampler->SetInterpolator ( nn_interpolator ); + //resampler->SetNumberOfThreads(1); + resampler->SetInput( itkImage ); resampler->SetTransform( base_transform ); resampler->SetReferenceImage( reference_image ); resampler->UseReferenceImageOn(); resampler->Update(); mitk::GrabItkImageMemory( resampler->GetOutput(), outputImage); } }; } // end namespace #endif // MITKPYRAMIDIMAGEREGISTRATION_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h index 7983bf02ea..8a772722fe 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h @@ -1,126 +1,160 @@ /*=================================================================== 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 MITKPYRAMIDREGISTRATIONMETHODHELPER_H #define MITKPYRAMIDREGISTRATIONMETHODHELPER_H #include #include #include -#include -#include +#include +#include +#include + +#include #include #include +#include +#include + #include #include #include "mitkImageAccessByItk.h" /** * @brief Provides same functionality as \a AccessTwoImagesFixedDimensionByItk for a subset of types * * For now, the subset defined is only short and float. */ #define AccessTwoImagesFixedDimensionTypeSubsetByItk(mitkImage1, mitkImage2, itkImageTypeFunction, dimension) \ { \ const mitk::PixelType& pixelType1 = mitkImage1->GetPixelType(); \ const mitk::PixelType& pixelType2 = mitkImage2->GetPixelType(); \ const mitk::Image* constImage1 = mitkImage1; \ const mitk::Image* constImage2 = mitkImage2; \ mitk::Image* nonConstImage1 = const_cast(constImage1); \ mitk::Image* nonConstImage2 = const_cast(constImage2); \ nonConstImage1->Update(); \ nonConstImage2->Update(); \ _checkSpecificDimension(mitkImage1, (dimension)); \ _checkSpecificDimension(mitkImage2, (dimension)); \ _accessTwoImagesByItkForEach(itkImageTypeFunction, ((short, dimension))((unsigned short, dimension))((float, dimension))((double, dimension)), ((short, dimension))((unsigned short, dimension))((float, dimension))((double, dimension)) ) \ { \ std::string msg("Pixel type "); \ msg.append(pixelType1.GetComponentTypeAsString() ); \ msg.append(" or pixel type "); \ msg.append(pixelType2.GetComponentTypeAsString() ); \ msg.append(" is not in " MITK_PP_STRINGIZE(MITK_ACCESSBYITK_TYPES_DIMN_SEQ(dimension))); \ throw mitk::AccessByItkException(msg); \ } \ } /** * @brief The PyramidOptControlCommand class stears the step lenght of the * gradient descent optimizer in multi-scale registration methods */ template class PyramidOptControlCommand : public itk::Command { public: typedef itk::RegularStepGradientDescentOptimizer OptimizerType; mitkClassMacro(PyramidOptControlCommand, itk::Command) itkFactorylessNewMacro(Self) itkCloneMacro(Self) void Execute(itk::Object *caller, const itk::EventObject & /*event*/) { RegistrationType* registration = dynamic_cast< RegistrationType* >( caller ); MITK_DEBUG << "\t - Pyramid level " << registration->GetCurrentLevel(); if( registration->GetCurrentLevel() == 0 ) return; OptimizerType* optimizer = dynamic_cast< OptimizerType* >(registration->GetOptimizer()); - MITK_INFO << optimizer->GetStopConditionDescription() << "\n" + MITK_INFO /*<< optimizer->GetStopConditionDescription() << "\n"*/ << optimizer->GetValue() << " : " << optimizer->GetCurrentPosition(); optimizer->SetMaximumStepLength( optimizer->GetMaximumStepLength() * 0.25f ); optimizer->SetMinimumStepLength( optimizer->GetMinimumStepLength() * 0.1f ); // optimizer->SetNumberOfIterations( optimizer->GetNumberOfIterations() * 1.5f ); } void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/){} }; template class OptimizerIterationCommand : public itk::Command { public: mitkClassMacro(OptimizerIterationCommand, itk::Command) itkFactorylessNewMacro(Self) itkCloneMacro(Self) void Execute(itk::Object *caller, const itk::EventObject & /*event*/) { + OptimizerType* optimizer = dynamic_cast< OptimizerType* >( caller ); unsigned int currentIter = optimizer->GetCurrentIteration(); MITK_INFO << "[" << currentIter << "] : " << optimizer->GetValue() << " : " << optimizer->GetCurrentPosition(); } void Execute(const itk::Object * /*object*/, const itk::EventObject & /*event*/) { } }; +template +class OptimizerIterationCommandv4 : public itk::Command +{ +public: + itkNewMacro( OptimizerIterationCommandv4 ) + + void Execute(itk::Object *object, const itk::EventObject & event) + { + OptimizerType* optimizer = dynamic_cast< OptimizerType* >( object ); + + if( typeid( event ) != typeid( itk::IterationEvent ) ) + { return; } + + unsigned int currentIter = optimizer->GetCurrentIteration(); + MITK_INFO << "[" << currentIter << "] : " << optimizer->GetCurrentMetricValue() << " : " + << optimizer->GetMetric()->GetParameters() ; + //<< " : " << optimizer->GetScales(); + + } + + void Execute(const itk::Object * object, const itk::EventObject & event) + { + + + } +}; + #endif // MITKPYRAMIDREGISTRATIONMETHODHELPER_H diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h rename to Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx similarity index 100% rename from Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx rename to Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h index 90c625658f..66ad3f8f24 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h @@ -1,146 +1,153 @@ /*=================================================================== 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 ITKSPLITDWIMAGEFILTER_H #define ITKSPLITDWIMAGEFILTER_H #include "itkImageToImageFilter.h" #include namespace itk { /** \class SplitDWImageFilter * * \brief Splits a DW-Image passed in as input into a 3D-t image where each volume coresponds to a * gradient image ( or the unweighted b0 ) * * Several applications require to get the gradient images as a separate volume, f.e. the registration for * head-motion correction. Also a reduction of the DW Image is possible when combined with its counterpart filter, * the \sa mitkImageToDiffusionImageSource, which can reinterpret a 3d+t (scalar) image into a diffusion weighted image. */ template< class TInputImagePixelType, class TOutputImagePixelType > class SplitDWImageFilter : public ImageToImageFilter< VectorImage, Image< TOutputImagePixelType, 4 > > { public: typedef SplitDWImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef TInputImagePixelType InputPixelType; typedef TOutputImagePixelType OutputPixelType; typedef ImageToImageFilter< VectorImage, Image< TOutputImagePixelType, 4 > > Superclass; /** typedefs from superclass */ typedef typename Superclass::InputImageType InputImageType; //typedef typename Superclass::OutputImageType OutputImageType; typedef Image< TOutputImagePixelType, 4 > OutputImageType; typedef typename OutputImageType::RegionType OutputImageRegionType; typedef std::vector< unsigned int > IndexListType; typedef std::map< unsigned int, IndexListType> BValueMapType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(SplitDWImageFilter, ImageToImageFilter); /** * @brief Set the indices of the images to be extracted. * * */ void SetExtractIndices( IndexListType list) { m_IndexList = list; } /** * @brief Extract all images * * The same as setting SetExtractIndices( ) with [0,1,2,...,N] */ void SetExtractAll() { m_ExtractAllImages = true; } /** * @brief Selects only the weighted images with b-value above the given b_threshold to be extracted * * Setting b_threshold to 0 will do the same as \sa SetExtractAll. Please note that some images have no true * unweighted images as the minimal b-value is something like 5 so for extracting all. * * @note It will reorder the images! * * @param b_threshold the minimal b-value to be extracted * @param map the map with b-values to the corresponding image * * @sa GetIndexList */ void SetExtractAllAboveThreshold( double b_threshold, BValueMapType map); + /** + * @brief SetExtractSingleShell + * @param b_value b-value of the shell to be extracted + * @param tol the tolerance of the shell choice, i.e. all shells within [ b_value - tol, b_value + tol ] will be extracted + */ + void SetExtractSingleShell( double b_value, BValueMapType map, double tol ); + /** * @brief Returns the index list used for extraction * * The list is necessary for further processing, especially when a b-value threshold is used ( like in \sa SetExtractAllAboveThreshold ) * * @return The index list used during the extraction */ const IndexListType GetIndexList() const { return m_IndexList; } protected: SplitDWImageFilter(); virtual ~SplitDWImageFilter(){}; void GenerateData(); /** The dimension of the output does not match the dimension of the input hence we need to re-implement the CopyInformation method to avoid executing the default implementation which tries to copy the input information to the output */ virtual void CopyInformation( const DataObject *data); /** Override of the ProcessObject::GenerateOutputInformation() because of different dimensionality of the input and the output */ virtual void GenerateOutputInformation(); IndexListType m_IndexList; bool m_ExtractAllImages; }; } //end namespace itk #ifndef ITK_MANUAL_INSTANTIATION #include "itkSplitDWImageFilter.txx" #endif #endif // ITKSPLITDWIMAGEFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx index 79a97ab912..338f2c5e4e 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.txx @@ -1,189 +1,212 @@ /*=================================================================== 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 __itkSplitDWImageFilter_txx #define __itkSplitDWImageFilter_txx #include "itkSplitDWImageFilter.h" #include #include template< class TInputImagePixelType, class TOutputImagePixelType> itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > ::SplitDWImageFilter() : m_IndexList(0), m_ExtractAllImages(true) { this->SetNumberOfRequiredInputs(1); } template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > ::CopyInformation( const DataObject* /*data*/) { } template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType >::GenerateOutputInformation() { } template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > ::SetExtractAllAboveThreshold(double b_threshold, BValueMapType map) { m_ExtractAllImages = false; m_IndexList.clear(); // create the index list following the given threshold // iterate over the b-value map BValueMapType::const_iterator bvalueIt = map.begin(); while( bvalueIt != map.end() ) { // check threshold if (bvalueIt->first > b_threshold) { // the map contains an index container, this needs to be inserted into the // index list IndexListType::const_iterator listIt = bvalueIt->second.begin(); while( listIt != bvalueIt->second.end() ) { m_IndexList.push_back( *listIt ); ++listIt; } } ++bvalueIt; } } +template< class TInputImagePixelType, + class TOutputImagePixelType> +void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType > +::SetExtractSingleShell(double b_value, BValueMapType map, double tol) +{ + m_ExtractAllImages = false; + m_IndexList.clear(); + + // create index list + BValueMapType::const_iterator bvalueIt = map.begin(); + while( bvalueIt != map.end() ) + { + IndexListType::const_iterator listIt = bvalueIt->second.begin(); + if( std::fabs( bvalueIt->first - b_value) < tol) + { + m_IndexList.insert( m_IndexList.begin(), bvalueIt->second.begin(), bvalueIt->second.end() ); + ++listIt; + } + + ++bvalueIt; + } +} + template< class TInputImagePixelType, class TOutputImagePixelType> void itk::SplitDWImageFilter< TInputImagePixelType, TOutputImagePixelType >::GenerateData() { if( m_IndexList.empty() && !m_ExtractAllImages ) { itkExceptionMacro(<<"The index list is empty and the option to extract all images is disabled. "); } // construct the index list (for each component) if( m_ExtractAllImages ) { m_IndexList.clear(); for( unsigned int i=0; i< this->GetInput()->GetNumberOfComponentsPerPixel(); i++) m_IndexList.push_back(i); } // declare the output image // this will have the b0 images stored as timesteps typename OutputImageType::Pointer outputImage = this->GetOutput(); //OutputImageType::New(); // allocate image with // - dimension: [DimX, DimY, DimZ, size(IndexList) ] // - spacing: old one, 1.0 time typename OutputImageType::SpacingType spacing; spacing.Fill(1); typename OutputImageType::PointType origin; origin.Fill(0); OutputImageRegionType outputRegion; // the spacing and origin corresponds to the respective values in the input image (3D) // the same for the region for (unsigned int i=0; i< 3; i++) { spacing[i] = (this->GetInput()->GetSpacing())[i]; origin[i] = (this->GetInput()->GetOrigin())[i]; outputRegion.SetSize(i, this->GetInput()->GetLargestPossibleRegion().GetSize()[i]); outputRegion.SetIndex(i, this->GetInput()->GetLargestPossibleRegion().GetIndex()[i]); } // number of timesteps = number of b0 images outputRegion.SetSize(3, m_IndexList.size()); outputRegion.SetIndex( 3, 0 ); // output image direction (4x4) typename OutputImageType::DirectionType outputDirection; //initialize to identity outputDirection.SetIdentity(); // get the input image direction ( 3x3 matrix ) typename InputImageType::DirectionType inputDirection = this->GetInput()->GetDirection(); for( unsigned int i=0; i< 3; i++) { outputDirection(0,i) = inputDirection(0,i); outputDirection(1,i) = inputDirection(1,i); outputDirection(2,i) = inputDirection(2,i); } outputImage->SetSpacing( spacing ); outputImage->SetOrigin( origin ); outputImage->SetDirection( outputDirection ); outputImage->SetRegions( outputRegion ); outputImage->Allocate(); // input iterator itk::ImageRegionConstIterator inputIt( this->GetInput(), this->GetInput()->GetLargestPossibleRegion() ); // we want to iterate separately over each 3D volume of the output image // so reset the regions last dimension outputRegion.SetSize(3,1); unsigned int currentTimeStep = 0; // for each index in the iterator list, extract the image and insert it as a new time step for(IndexListType::const_iterator indexIt = m_IndexList.begin(); indexIt != m_IndexList.end(); indexIt++) { // set the time step outputRegion.SetIndex(3, currentTimeStep); itk::ImageRegionIterator< OutputImageType> outputIt( outputImage.GetPointer(), outputRegion ); // iterate over the current b0 image and store it to corresponding place outputIt.GoToBegin(); inputIt.GoToBegin(); while( !outputIt.IsAtEnd() && !inputIt.IsAtEnd() ) { // the input vector typename InputImageType::PixelType vec = inputIt.Get(); outputIt.Set( vec[*indexIt]); ++outputIt; ++inputIt; } // increase time step currentTimeStep++; } } #endif // __itkSplitDWImageFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h index a635a202dd..cbed2d792b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h @@ -1,158 +1,168 @@ /*=================================================================== 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_TensorImageToDiffusionImageFilter_h_ #define _itk_TensorImageToDiffusionImageFilter_h_ #include "itkImageToImageFilter.h" -#include +#include #include namespace itk { template class TensorImageToDiffusionImageFilter : public ImageToImageFilter,3>, itk::VectorImage > { public: typedef TInputScalarType InputScalarType; typedef itk::DiffusionTensor3D InputPixelType; typedef itk::Image InputImageType; typedef typename InputImageType::RegionType InputImageRegionType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef typename OutputImageType::RegionType OutputImageRegionType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef typename BaselineImageType::RegionType BaselineImageRegionType; + typedef itk::Image< short, 3> MaskImageType; + typedef MaskImageType::RegionType MaskImageRegionType; + typedef TensorImageToDiffusionImageFilter Self; typedef ImageToImageFilter Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; itkTypeMacro (TensorImageToDiffusionImageFilter, ImageToImageFilter); itkStaticConstMacro (ImageDimension, unsigned int, OutputImageType::ImageDimension); itkFactorylessNewMacro(Self) itkCloneMacro(Self) - typedef vnl_vector_fixed< double, 3 > GradientType; - typedef itk::VectorContainer< unsigned int, GradientType > GradientListType; + typedef vnl_vector_fixed GradientType; + typedef VectorContainer GradientListType; + typedef GradientListType::Pointer GradientListPointerType; /** Manually Set/Get a list of gradients */ - void SetGradientList(GradientListType::Pointer list) + void SetGradientList(const GradientListPointerType list) { m_GradientList = list; this->Modified(); } - GradientListType GetGradientList(void) const + GradientListPointerType GetGradientList(void) const {return m_GradientList;} void SetBValue( const double& bval) { m_BValue = bval; } + void SetMaskImage( MaskImageType::Pointer maskimage ) + { + m_MaskImage = maskimage; + } + /** * @brief Set an external baseline image for signal generation (optional) * * An option to enforce a specific baseline image. If none provided (default) the filter uses * the itk::TensorToL2NormImageFilter to generate the modelled baseline image. */ void SetExternalBaselineImage( typename BaselineImageType::Pointer bimage) { m_BaselineImage = bimage; } itkSetMacro(Min, OutputScalarType); itkSetMacro(Max, OutputScalarType); protected: TensorImageToDiffusionImageFilter() { m_BValue = 1.0; m_BaselineImage = 0; m_Min = 0.0; m_Max = 10000.0; } - ~TensorImageToDiffusionImageFilter(){} + virtual ~TensorImageToDiffusionImageFilter(){} void PrintSelf (std::ostream& os, Indent indent) const { Superclass::PrintSelf (os, indent); } - void BeforeThreadedGenerateData( void ); + virtual void BeforeThreadedGenerateData( void ); - void ThreadedGenerateData( const + virtual void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); //void GenerateData(); - - private: - - TensorImageToDiffusionImageFilter (const Self&); - void operator=(const Self&); - - GradientListType::Pointer m_GradientList; + GradientListPointerType m_GradientList; double m_BValue; typename BaselineImageType::Pointer m_BaselineImage; OutputScalarType m_Min; OutputScalarType m_Max; + MaskImageType::Pointer m_MaskImage; + + private: + + TensorImageToDiffusionImageFilter (const Self&); + void operator=(const Self&); + }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkTensorImageToDiffusionImageFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx index 4029ca4a0b..e812dfb51a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx @@ -1,211 +1,235 @@ /*=================================================================== 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. ===================================================================*/ /*========================================================================= 2 3 Program: Tensor ToolKit - TTK 4 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $ 5 Language: C++ 6 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ 7 Version: $Revision: 68 $ 8 9 Copyright (c) INRIA 2010. All rights reserved. 10 See LICENSE.txt for details. 11 12 This software is distributed WITHOUT ANY WARRANTY; without even 13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 14 PURPOSE. See the above copyright notices for more information. 15 16 =========================================================================*/ #ifndef _itk_TensorImageToDiffusionImageFilter_txx_ #define _itk_TensorImageToDiffusionImageFilter_txx_ #endif #include "itkTensorImageToDiffusionImageFilter.h" #include "itkTensorToL2NormImageFilter.h" #include "itkRescaleIntensityImageFilter.h" #include #include namespace itk { template void TensorImageToDiffusionImageFilter ::BeforeThreadedGenerateData() { if( m_GradientList->Size()==0 ) { throw itk::ExceptionObject (__FILE__,__LINE__,"Error: gradient list is empty, cannot generate DWI."); } if( m_BaselineImage.IsNull() ) { // create a B0 image by taking the norm of the tensor field * scale: typedef itk::TensorToL2NormImageFilter > TensorToL2NormFilterType; typename TensorToL2NormFilterType::Pointer myFilter1 = TensorToL2NormFilterType::New(); myFilter1->SetInput (this->GetInput()); try { myFilter1->Update(); } catch (itk::ExceptionObject &e) { std::cerr << e; return; } typename itk::RescaleIntensityImageFilter< itk::Image, BaselineImageType>::Pointer rescaler= itk::RescaleIntensityImageFilter, BaselineImageType>::New(); rescaler->SetOutputMinimum ( m_Min ); rescaler->SetOutputMaximum ( m_Max ); rescaler->SetInput ( myFilter1->GetOutput() ); try { rescaler->Update(); } catch (itk::ExceptionObject &e) { std::cerr << e; return; } m_BaselineImage = rescaler->GetOutput(); } typename OutputImageType::Pointer outImage = OutputImageType::New(); outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion()); outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetVectorLength(m_GradientList->Size()); outImage->Allocate(); this->SetNumberOfRequiredOutputs (1); this->SetNthOutput (0, outImage); } template -void -TensorImageToDiffusionImageFilter +void TensorImageToDiffusionImageFilter ::ThreadedGenerateData (const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId ) { typedef ImageRegionIterator IteratorOutputType; typedef ImageRegionConstIterator IteratorInputType; typedef ImageRegionConstIterator IteratorBaselineType; unsigned long numPixels = outputRegionForThread.GetNumberOfPixels(); unsigned long step = numPixels/100; unsigned long progress = 0; IteratorOutputType itOut (this->GetOutput(), outputRegionForThread); IteratorInputType itIn (this->GetInput(), outputRegionForThread); IteratorBaselineType itB0 (m_BaselineImage, outputRegionForThread); + typedef ImageRegionConstIterator< MaskImageType > IteratorMaskImageType; + IteratorMaskImageType itMask; + + if( m_MaskImage.IsNotNull() ) + { + itMask = IteratorMaskImageType( m_MaskImage, outputRegionForThread); + itMask.GoToBegin(); + } + if( threadId==0 ) { this->UpdateProgress (0.0); } - while(!itIn.IsAtEnd()) { if( this->GetAbortGenerateData() ) { throw itk::ProcessAborted(__FILE__,__LINE__); } InputPixelType T = itIn.Get(); BaselinePixelType b0 = itB0.Get(); OutputPixelType out; out.SetSize(m_GradientList->Size()); out.Fill(0); + short maskvalue = 1; + if( m_MaskImage.IsNotNull() ) + { + maskvalue = itMask.Get(); + ++itMask; + } + + std::vector b0_indices; if( b0 > 0) { - for( unsigned int i=0; iSize()-1; i++) + for( unsigned int i=0; iSize(); i++) { - - GradientType g = m_GradientList->at(i+1); + GradientType g = m_GradientList->at(i); // normalize vector so the following computations work const double twonorm = g.two_norm(); + if( twonorm < vnl_math::eps ) + { + b0_indices.push_back(i); + continue; + } + GradientType gn = g.normalize(); InputPixelType S; S[0] = gn[0]*gn[0]; S[1] = gn[1]*gn[0]; S[2] = gn[2]*gn[0]; S[3] = gn[1]*gn[1]; S[4] = gn[2]*gn[1]; S[5] = gn[2]*gn[2]; const double res = T[0]*S[0] + 2 * T[1]*S[1] + T[3]*S[3] + 2 * T[2]*S[2] + 2 * T[4]*S[4] + T[5]*S[5]; // check for corrupted tensor if (res>=0) { // estimate the bvalue from the base value and the norm of the gradient // - because of this estimation the vector have to be normalized beforehand // otherwise the modelled signal is wrong ( i.e. not scaled properly ) - const double bval = m_BValue * twonorm; - out[i+1] = static_cast( 1.0 * b0 * exp ( -1.0 * bval * res ) ); + const double bval = m_BValue * twonorm * twonorm; + out[i] = static_cast( maskvalue * 1.0 * b0 * exp ( -1.0 * bval * res ) ); } } } - out[0] = b0; + for(unsigned int idx = 0; idx < b0_indices.size(); idx++ ) + { + out[b0_indices.at(idx)] = b0; + } + itOut.Set(out); if( threadId==0 && step>0) { if( (progress%step)==0 ) { this->UpdateProgress ( double(progress)/double(numPixels) ); } } ++progress; ++itB0; ++itIn; ++itOut; } if( threadId==0 ) { this->UpdateProgress (1.0); } } } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt index 35a29c518f..517a120410 100644 --- a/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt @@ -1,10 +1,11 @@ MITK_CREATE_MODULE( SUBPROJECTS MITK-DTI INCLUDE_DIRS Algorithms Algorithms/Reconstruction Algorithms/Registration Algorithms/Reconstruction/MultishellProcessing DicomImport IODataStructures/DiffusionWeightedImages IODataStructures/QBallImages IODataStructures/TensorImages IODataStructures Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS MitkMapperExt MitkPlanarFigure MitkImageExtraction MitkSceneSerializationBase PACKAGE_DEPENDS VTK|vtkFiltersProgrammable ITK|ITKDistanceMap+ITKRegistrationCommon+ITKLabelVoting+ITKVTK Boost + ITK|ITKMetricsv4+ITKRegistrationMethodsv4 WARNINGS_AS_ERRORS ) add_subdirectory(Testing) diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp index 6dcbfbba98..50ce16c379 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp @@ -1,143 +1,145 @@ /*=================================================================== 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 __mitk_ImageToDiffusionImageSource_txx #define __mitk_ImageToDiffusionImageSource_txx #include "mitkImageToDiffusionImageSource.h" #include #include #include "mitkImageTimeSelector.h" #include #include #include "mitkImageCast.h" #include "mitkException.h" template< typename TPixelType> mitk::ImageToDiffusionImageSource::ImageToDiffusionImageSource() : m_SourceImage(0), m_GradientDirections(0) { } template< typename TPixelType> void mitk::ImageToDiffusionImageSource ::GenerateOutputInformation() { // sanity checks ( input not null, timesteps corresponds to number of directions specified ) if( m_GradientDirections.IsNull() || m_GradientDirections->empty() ) { mitkThrow() << "No gradient directions were set. Cannot proceed."; } if( m_SourceImage.IsNull() ) { mitkThrow() << "No input image was set. Cannot proceed."; } if( m_GradientDirections->size() != m_SourceImage->GetTimeSteps() ) { mitkThrow() << "Size mismatch between container size " << m_GradientDirections->size() << "and image volumes."<< m_SourceImage->GetTimeSteps() <<" Cannot proceed."; } // already pass in the meta-data typename OutputType::Pointer metaImage = OutputType::New(); // set identity as measurement frame vnl_matrix_fixed< double, 3, 3 > measurement_frame; measurement_frame.set_identity(); metaImage->SetMeasurementFrame( measurement_frame ); // set directions and bvalue metaImage->SetDirections(this->m_GradientDirections); metaImage->SetReferenceBValue(this->m_BValue); m_OutputCache = metaImage; m_CacheTime.Modified(); } template< typename TPixelType> void mitk::ImageToDiffusionImageSource ::GenerateData() { if ( ( ! m_OutputCache ) || ( this->GetMTime( ) > m_CacheTime.GetMTime( ) ) ) { this->GenerateOutputInformation(); itkWarningMacro("Cache regenerated!"); } // now cast the mitk image to the vector image and pass in as volume typedef itk::Image< TPixelType, 4 > InputItkType; typedef itk::Image< TPixelType, 3> SingleVolumeType; typedef itk::VectorImage< TPixelType, 3> VectorImageType; typedef itk::ComposeImageFilter< SingleVolumeType > ComposeFilterType; typename ComposeFilterType::Pointer vec_composer = ComposeFilterType::New(); mitk::ImageTimeSelector::Pointer t_selector = mitk::ImageTimeSelector::New(); t_selector->SetInput( m_SourceImage ); for( unsigned int i=0; i< m_SourceImage->GetTimeSteps(); i++) { t_selector->SetTimeNr(i); t_selector->Update(); typename SingleVolumeType::Pointer singleImageItk; mitk::CastToItkImage( t_selector->GetOutput(), singleImageItk ); vec_composer->SetInput( i, singleImageItk ); } try { vec_composer->Update(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Caugt exception while updating compose filter: " << e.what(); } m_OutputCache->SetVectorImage( vec_composer->GetOutput() ); // transfer to the output image static_cast(this->GetOutput()) ->SetVectorImage(m_OutputCache->GetVectorImage()); static_cast(this->GetOutput()) ->SetReferenceBValue(m_OutputCache->GetReferenceBValue()); static_cast(this->GetOutput()) ->SetMeasurementFrame(m_OutputCache->GetMeasurementFrame()); static_cast(this->GetOutput()) ->SetDirections(m_OutputCache->GetDirections()); + static_cast(this->GetOutput()) + ->UpdateBValueMap(); static_cast(this->GetOutput()) ->InitializeFromVectorImage(); } #endif //__mitk_ImageToDiffusionImageSource_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake b/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake index e1e74266f2..4b3b651899 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake @@ -1,12 +1,14 @@ set(MODULE_TESTS mitkFactoryRegistrationTest.cpp mitkDiffusionImageEqualTest.cpp mitkNonLocalMeansDenoisingTest.cpp ) set(MODULE_CUSTOM_TESTS mitkPyramidImageRegistrationMethodTest.cpp mitkDWHeadMotionCorrectionTest.cpp mitkImageReconstructionTest.cpp + mitkConvertDWITypeTest.cpp + mitkExtractSingleShellTest.cpp ) diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp similarity index 77% copy from Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp copy to Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp index ad3c899314..6998c8a3b3 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp @@ -1,57 +1,59 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "mitkIOUtil.h" #include "mitkDWIHeadMotionCorrectionFilter.h" +#include "mitkNrrdDiffusionImageWriter.h" typedef short DiffusionPixelType; typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; /** * @brief Custom test to provide CMD-line access to the mitk::DWIHeadMotionCorrectionFilter * * @param argv : Input and Output image full path */ -int mitkDWHeadMotionCorrectionTest( int argc, char* argv[] ) +int mitkConvertDWITypeTest( int argc, char* argv[] ) { - MITK_TEST_BEGIN("mitkDWHeadMotionCorrectionTest"); + MITK_TEST_BEGIN("mitkConvertDWITypeTest"); MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output."); mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); DiffusionImageType* dwimage = static_cast( inputImage.GetPointer() ); - mitk::DWIHeadMotionCorrectionFilter::Pointer corrfilter = - mitk::DWIHeadMotionCorrectionFilter::New(); - corrfilter->SetInput( dwimage ); - corrfilter->Update(); + mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::Pointer dwiwriter = + mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::New(); + + dwiwriter->SetInput( dwimage ); + dwiwriter->SetFileName( argv[2] ); try { - mitk::IOUtil::SaveBaseData(corrfilter->GetOutput(), argv[2]); + dwiwriter->Update(); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Catched exception: " << e.what(); mitkThrow() << "Failed with exception from subprocess!"; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp index ad3c899314..3771fe793b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp @@ -1,57 +1,60 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "mitkIOUtil.h" #include "mitkDWIHeadMotionCorrectionFilter.h" typedef short DiffusionPixelType; typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; /** * @brief Custom test to provide CMD-line access to the mitk::DWIHeadMotionCorrectionFilter * * @param argv : Input and Output image full path */ int mitkDWHeadMotionCorrectionTest( int argc, char* argv[] ) { MITK_TEST_BEGIN("mitkDWHeadMotionCorrectionTest"); MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output."); + itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1); + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); DiffusionImageType* dwimage = static_cast( inputImage.GetPointer() ); mitk::DWIHeadMotionCorrectionFilter::Pointer corrfilter = mitk::DWIHeadMotionCorrectionFilter::New(); corrfilter->SetInput( dwimage ); + corrfilter->SetAverageUnweighted(false); corrfilter->Update(); try { mitk::IOUtil::SaveBaseData(corrfilter->GetOutput(), argv[2]); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Catched exception: " << e.what(); mitkThrow() << "Failed with exception from subprocess!"; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp new file mode 100644 index 0000000000..93aa88088e --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp @@ -0,0 +1,103 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkTestingMacros.h" +#include "mitkIOUtil.h" + +#include + +#include "mitkDWIHeadMotionCorrectionFilter.h" +#include "mitkNrrdDiffusionImageWriter.h" + +typedef short DiffusionPixelType; +typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; + + +int mitkExtractSingleShellTest( int argc, char* argv[] ) +{ + MITK_TEST_BEGIN("mitkExtractSingleShellTest"); + + MITK_TEST_CONDITION_REQUIRED( argc > 3, "Specify input and output and the shell to be extracted"); + + /* + 1. Get input data + */ + mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); + DiffusionImageType* dwimage = + static_cast( inputImage.GetPointer() ); + + MITK_TEST_CONDITION_REQUIRED( dwimage != NULL, "Input is a dw-image"); + + unsigned int extract_value = 0; + std::istringstream input(argv[3]); + + input >> extract_value; + + typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter FilterType; + typedef DiffusionImageType::BValueMap BValueMap; + + // GetShellSelection from GUI + BValueMap shellSelectionMap; + BValueMap originalShellMap = dwimage->GetB_ValueMap(); + std::vector newNumGradientDirections; + + shellSelectionMap[extract_value] = originalShellMap[extract_value]; + newNumGradientDirections.push_back( originalShellMap[extract_value].size() ) ; + + DiffusionImageType::GradientDirectionContainerType::Pointer gradientContainer = dwimage->GetDirections(); + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(dwimage->GetVectorImage()); + filter->SetOriginalGradientDirections(gradientContainer); + filter->SetNumGradientDirections(newNumGradientDirections); + filter->SetOriginalBValueMap(originalShellMap); + filter->SetShellSelectionBValueMap(shellSelectionMap); + + try + { + filter->Update(); + } + catch( const itk::ExceptionObject& e) + { + MITK_TEST_FAILED_MSG( << "Failed due to ITK exception: " << e.what() ); + } + + DiffusionImageType::Pointer image = DiffusionImageType::New(); + image->SetVectorImage( filter->GetOutput() ); + image->SetB_Value(dwimage->GetB_Value()); + image->SetDirections(filter->GetGradientDirections()); + image->SetMeasurementFrame(dwimage->GetMeasurementFrame()); + image->InitializeFromVectorImage(); + /* + * 3. Write output data + **/ + mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::Pointer dwiwriter = + mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::New(); + + dwiwriter->SetInput( image ); + dwiwriter->SetFileName( argv[2] ); + + try + { + dwiwriter->Update(); + } + catch( const itk::ExceptionObject& e) + { + MITK_ERROR << "Catched exception: " << e.what(); + mitkThrow() << "Failed with exception from subprocess!"; + } + + MITK_TEST_END(); +} diff --git a/Modules/DiffusionImaging/DiffusionCore/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake index ac90de0f2d..5ff5e61b19 100644 --- a/Modules/DiffusionImaging/DiffusionCore/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake @@ -1,129 +1,131 @@ set(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp # DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkQBallImage.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImage.cpp #IODataStructures/mitkRegistrationObject.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/itkDwiGradientLengthCorrectionFilter.cpp + Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h # Registration Algorithms & Co. Algorithms/Registration/mitkRegistrationWrapper.cpp Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp + # Algorithms/Registration/mitkRegistrationMethodITK4.cpp # MultishellProcessing Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkOdfVtkMapper2D.h # Reconstruction Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Algorithms/Reconstruction/itkPointShell.h Algorithms/Reconstruction/itkOrientationDistributionFunction.h Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h # MultishellProcessing Algorithms/Reconstruction/MultishellProcessing/itkRadialMultishellToSingleshellImageFilter.h Algorithms/Reconstruction/MultishellProcessing/itkDWIVoxelFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h Algorithms/itkMergeDiffusionImagesFilter.h Algorithms/itkDwiPhantomGenerationFilter.h Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h Algorithms/itkMrtrixPeakImageConverter.h Algorithms/itkFslPeakImageConverter.h Algorithms/itkShCoefficientImageImporter.h Algorithms/itkShCoefficientImageExporter.h Algorithms/itkOdfMaximaExtractionFilter.h Algorithms/itkResampleDwiImageFilter.h Algorithms/itkDwiGradientLengthCorrectionFilter.h Algorithms/itkAdcImageFilter.h Algorithms/itkSplitDWImageFilter.h Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h Algorithms/mitkDiffusionImageToDiffusionImageFilter.h Algorithms/itkNonLocalMeansDenoisingFilter.h Algorithms/itkVectorImageToImageFilter.h ) set( TOOL_FILES ) diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 05776f0ff2..0c89c03aaf 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,78 +1,79 @@ set(CPP_FILES ## IO datastructures IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkTrackvis.cpp IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Tractography Algorithms/GibbsTracking/mitkParticleGrid.cpp Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp Algorithms/GibbsTracking/mitkEnergyComputer.cpp Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp Algorithms/GibbsTracking/mitkFiberBuilder.cpp Algorithms/GibbsTracking/mitkSphereInterpolator.cpp ) set(H_FILES # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkTrackvis.h IODataStructures/mitkFiberfoxParameters.h # Algorithms Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h - Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h + # moved to DiffusionCore + #Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h Algorithms/itkFibersFromPlanarFiguresFilter.h Algorithms/itkTractsToDWIImageFilter.h Algorithms/itkTractsToVectorImageFilter.h Algorithms/itkKspaceImageFilter.h Algorithms/itkDftImageFilter.h Algorithms/itkAddArtifactsToDwiImageFilter.h Algorithms/itkFieldmapGeneratorFilter.h Algorithms/itkEvaluateDirectionImagesFilter.h Algorithms/itkEvaluateTractogramDirectionsFilter.h # (old) Tractography Algorithms/itkGibbsTrackingFilter.h Algorithms/itkStochasticTractographyFilter.h Algorithms/itkStreamlineTrackingFilter.h Algorithms/GibbsTracking/mitkParticle.h Algorithms/GibbsTracking/mitkParticleGrid.h Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.h Algorithms/GibbsTracking/mitkSimpSamp.h Algorithms/GibbsTracking/mitkEnergyComputer.h Algorithms/GibbsTracking/mitkGibbsEnergyComputer.h Algorithms/GibbsTracking/mitkSphereInterpolator.h Algorithms/GibbsTracking/mitkFiberBuilder.h # Signal Models SignalModels/mitkDiffusionSignalModel.h SignalModels/mitkTensorModel.h SignalModels/mitkBallModel.h SignalModels/mitkDotModel.h SignalModels/mitkAstroStickModel.h SignalModels/mitkStickModel.h SignalModels/mitkDiffusionNoiseModel.h SignalModels/mitkRicianNoiseModel.h SignalModels/mitkChiSquareNoiseModel.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin # Shaders Shaders/mitkShaderFiberClipping.xml ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp index 1df5dff8ef..6df9416faa 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp @@ -1,1022 +1,1025 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkTensorReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include #include #include #include #include // itk includes #include "itkTimeProbe.h" //#include "itkTensor.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "mitkTeemDiffusionTensor3DReconstructionImageFilter.h" #include "itkDiffusionTensor3DReconstructionImageFilter.h" #include "itkTensorImageToDiffusionImageFilter.h" #include "itkPointShell.h" #include "itkVector.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkTensorReconstructionWithEigenvalueCorrectionFilter.h" //#include "itkFreeWaterEliminationFilter.h" #include "mitkProperties.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "mitkDiffusionImageMapper.h" #include "mitkLookupTableProperty.h" #include "mitkLookupTable.h" #include "mitkImageStatisticsHolder.h" #include #include #include #include const std::string QmitkTensorReconstructionView::VIEW_ID = "org.mitk.views.tensorreconstruction"; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; using namespace berry; QmitkTensorReconstructionView::QmitkTensorReconstructionView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL) { m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); } QmitkTensorReconstructionView::~QmitkTensorReconstructionView() { } void QmitkTensorReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkTensorReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); Advanced1CheckboxClicked(); } } void QmitkTensorReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkTensorReconstructionView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkTensorReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_StartReconstruction), SIGNAL(clicked()), this, SLOT(Reconstruct()) ); connect( (QObject*)(m_Controls->m_Advanced1), SIGNAL(clicked()), this, SLOT(Advanced1CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_TensorsToDWIButton), SIGNAL(clicked()), this, SLOT(TensorsToDWI()) ); connect( (QObject*)(m_Controls->m_TensorsToQbiButton), SIGNAL(clicked()), this, SLOT(TensorsToQbi()) ); connect( (QObject*)(m_Controls->m_ResidualButton), SIGNAL(clicked()), this, SLOT(ResidualCalculation()) ); connect( (QObject*)(m_Controls->m_PerSliceView), SIGNAL(pointSelected(int, int)), this, SLOT(ResidualClicked(int, int)) ); } } void QmitkTensorReconstructionView::ResidualClicked(int slice, int volume) { // Use image coord to reset crosshair // Find currently selected diffusion image // Update Label // to do: This position should be modified in order to skip B0 volumes that are not taken into account // when calculating residuals // Find the diffusion image mitk::DiffusionImage* diffImage; mitk::DataNode::Pointer correctNode; mitk::Geometry3D* geometry; if (m_DiffusionImage.IsNotNull()) { diffImage = static_cast*>(m_DiffusionImage->GetData()); geometry = diffImage->GetGeometry(); // Remember the node whose display index must be updated correctNode = mitk::DataNode::New(); correctNode = m_DiffusionImage; } if(diffImage != NULL) { typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; GradientDirectionContainerType::Pointer dirs = diffImage->GetDirections(); for(unsigned int i=0; iSize() && i<=volume; i++) { GradientDirectionType grad = dirs->ElementAt(i); // check if image is b0 weighted if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001) { volume++; } } QString pos = "Volume: "; pos.append(QString::number(volume)); pos.append(", Slice: "); pos.append(QString::number(slice)); m_Controls->m_PositionLabel->setText(pos); if(correctNode) { int oldDisplayVal; correctNode->GetIntProperty("DisplayChannel", oldDisplayVal); std::string oldVal = QString::number(oldDisplayVal).toStdString(); std::string newVal = QString::number(volume).toStdString(); correctNode->SetIntProperty("DisplayChannel",volume); correctNode->SetSelected(true); this->FirePropertyChanged("DisplayChannel", oldVal, newVal); correctNode->UpdateOutputInformation(); mitk::Point3D p3 = m_MultiWidget->GetCrossPosition(); itk::Index<3> ix; geometry->WorldToIndex(p3, ix); // ix[2] = slice; mitk::Vector3D vec; vec[0] = ix[0]; vec[1] = ix[1]; vec[2] = slice; mitk::Vector3D v3New; geometry->IndexToWorld(vec, v3New); mitk::Point3D origin = geometry->GetOrigin(); mitk::Point3D p3New; p3New[0] = v3New[0] + origin[0]; p3New[1] = v3New[1] + origin[1]; p3New[2] = v3New[2] + origin[2]; m_MultiWidget->MoveCrossToPosition(p3New); m_MultiWidget->RequestUpdate(); } } } void QmitkTensorReconstructionView::Advanced1CheckboxClicked() { bool check = m_Controls-> m_Advanced1->isChecked(); m_Controls->frame->setVisible(check); } void QmitkTensorReconstructionView::Activated() { QmitkFunctionality::Activated(); } void QmitkTensorReconstructionView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkTensorReconstructionView::ResidualCalculation() { // Extract dwi and dti from current selection // In case of multiple selections, take the first one, since taking all combinations is not meaningful mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); mitk::DiffusionImage::Pointer diffImage = mitk::DiffusionImage::New(); TensorImageType::Pointer tensorImage; std::string nodename; if(m_DiffusionImage.IsNotNull()) { diffImage = static_cast*>(m_DiffusionImage->GetData()); } else return; if(m_TensorImage.IsNotNull()) { mitk::TensorImage* mitkVol; mitkVol = static_cast(m_TensorImage->GetData()); mitk::CastToItkImage(mitkVol, tensorImage); m_TensorImage->GetStringProperty("name", nodename); } else return; typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; mitk::DiffusionImage::GradientDirectionContainerType* gradients = diffImage->GetDirections(); // Find the min and the max values from a baseline image mitk::ImageStatisticsHolder *stats = diffImage->GetStatistics(); //Initialize filter that calculates the modeled diffusion weighted signals FilterType::Pointer filter = FilterType::New(); filter->SetInput( tensorImage ); filter->SetBValue(diffImage->GetReferenceBValue()); filter->SetGradientList(gradients); filter->SetMin(stats->GetScalarValueMin()); filter->SetMax(stats->GetScalarValueMax()); filter->Update(); // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetReferenceBValue(diffImage->GetReferenceBValue()); image->SetDirections(gradients); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dwi"); node->SetName(newname.toAscii()); GetDefaultDataStorage()->Add(node); mitk::DiffusionImage::BValueMap map =image->GetBValueMap(); mitk::DiffusionImage::IndicesVector b0Indices = map[0]; typedef itk::ResidualImageFilter ResidualImageFilterType; ResidualImageFilterType::Pointer residualFilter = ResidualImageFilterType::New(); residualFilter->SetInput(diffImage->GetVectorImage()); residualFilter->SetSecondDiffusionImage(image->GetVectorImage()); residualFilter->SetGradients(gradients); residualFilter->SetB0Index(b0Indices[0]); residualFilter->SetB0Threshold(30); residualFilter->Update(); itk::Image::Pointer residualImage = itk::Image::New(); residualImage = residualFilter->GetOutput(); mitk::Image::Pointer mitkResImg = mitk::Image::New(); mitk::CastToMitkImage(residualImage, mitkResImg); stats = mitkResImg->GetStatistics(); float min = stats->GetScalarValueMin(); float max = stats->GetScalarValueMax(); mitk::LookupTableProperty::Pointer lutProp = mitk::LookupTableProperty::New(); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(min, max); // If you don't want to use the whole color range, you can use // SetValueRange, SetHueRange, and SetSaturationRange lookupTable->Build(); vtkSmartPointer reversedlookupTable = vtkSmartPointer::New(); reversedlookupTable->SetTableRange(min+1, max); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookupTable->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } lut->SetVtkLookupTable(lookupTable); lutProp->SetLookupTable(lut); // Create lookuptable mitk::DataNode::Pointer resNode=mitk::DataNode::New(); resNode->SetData( mitkResImg ); resNode->SetName("Residual Image"); resNode->SetProperty("LookupTable", lutProp); bool b; resNode->GetBoolProperty("use color", b); resNode->SetBoolProperty("use color", false); GetDefaultDataStorage()->Add(resNode); m_MultiWidget->RequestUpdate(); // Draw Graph std::vector means = residualFilter->GetMeans(); std::vector q1s = residualFilter->GetQ1(); std::vector q3s = residualFilter->GetQ3(); std::vector percentagesOfOUtliers = residualFilter->GetPercentagesOfOutliers(); m_Controls->m_ResidualAnalysis->SetMeans(means); m_Controls->m_ResidualAnalysis->SetQ1(q1s); m_Controls->m_ResidualAnalysis->SetQ3(q3s); m_Controls->m_ResidualAnalysis->SetPercentagesOfOutliers(percentagesOfOUtliers); if(m_Controls->m_PercentagesOfOutliers->isChecked()) { m_Controls->m_ResidualAnalysis->DrawPercentagesOfOutliers(); } else { m_Controls->m_ResidualAnalysis->DrawMeans(); } // Draw Graph for volumes per slice in the QGraphicsView std::vector< std::vector > outliersPerSlice = residualFilter->GetOutliersPerSlice(); int xSize = outliersPerSlice.size(); if(xSize == 0) { return; } int ySize = outliersPerSlice[0].size(); // Find maximum in outliersPerSlice double maxOutlier= 0.0; for(int i=0; imaxOutlier) { maxOutlier = outliersPerSlice[i][j]; } } } // Create some QImage QImage qImage(xSize, ySize, QImage::Format_RGB32); QImage legend(1, 256, QImage::Format_RGB32); QRgb value; vtkSmartPointer lookup = vtkSmartPointer::New(); lookup->SetTableRange(0.0, maxOutlier); lookup->Build(); reversedlookupTable->SetTableRange(0, maxOutlier); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookup->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } // Fill qImage for(int i=0; iMapValue(out); int r, g, b; r = _rgba[0]; g = _rgba[1]; b = _rgba[2]; value = qRgb(r, g, b); qImage.setPixel(i,j,value); } } for(int i=0; i<256; i++) { double* rgba = lookup->GetTableValue(i); int r, g, b; r = rgba[0]*255; g = rgba[1]*255; b = rgba[2]*255; value = qRgb(r, g, b); legend.setPixel(0,255-i,value); } QString upper = QString::number(maxOutlier, 'g', 3); upper.append(" %"); QString lower = QString::number(0.0); lower.append(" %"); m_Controls->m_UpperLabel->setText(upper); m_Controls->m_LowerLabel->setText(lower); QGraphicsScene* scene = new QGraphicsScene; QGraphicsScene* scene2 = new QGraphicsScene; QPixmap pixmap(QPixmap::fromImage(qImage)); QGraphicsPixmapItem *item = new QGraphicsPixmapItem( pixmap, 0, scene); item->scale(10.0, 3.0); QPixmap pixmap2(QPixmap::fromImage(legend)); QGraphicsPixmapItem *item2 = new QGraphicsPixmapItem( pixmap2, 0, scene2); item2->scale(20.0, 1.0); m_Controls->m_PerSliceView->SetResidualPixmapItem(item); m_Controls->m_PerSliceView->setScene(scene); m_Controls->m_LegendView->setScene(scene2); m_Controls->m_PerSliceView->show(); m_Controls->m_PerSliceView->repaint(); m_Controls->m_LegendView->setHorizontalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->setVerticalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->show(); m_Controls->m_LegendView->repaint(); } void QmitkTensorReconstructionView::Reconstruct() { int method = m_Controls->m_ReconctructionMethodBox->currentIndex(); switch (method) { case 0: ItkTensorReconstruction(m_DiffusionImages); break; case 1: TensorReconstructionWithCorr(m_DiffusionImages); break; default: ItkTensorReconstruction(m_DiffusionImages); } } void QmitkTensorReconstructionView::TensorReconstructionWithCorr (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientDirectionContainerType; DiffusionImageType* vols = static_cast((*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // TENSOR RECONSTRUCTION clock.Start(); MITK_INFO << "Tensor reconstruction with correction for negative eigenvalues"; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::TensorReconstructionWithEigenvalueCorrectionFilter< DiffusionPixelType, TTensorPixelType > ReconstructionFilter; float b0Threshold = m_Controls->m_TensorReconstructionThreshold->value(); GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = vols->GetDirections()->Begin(); it != vols->GetDirections()->End(); it++) { gradientContainerCopy->push_back(it.Value()); } ReconstructionFilter::Pointer reconFilter = ReconstructionFilter::New(); reconFilter->SetGradientImage( gradientContainerCopy, vols->GetVectorImage() ); reconFilter->SetBValue(vols->GetReferenceBValue()); reconFilter->SetB0Threshold(b0Threshold); reconFilter->Update(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer outputTensorImg = reconFilter->GetOutput(); typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(outputTensorImg, outputTensorImg->GetRequestedRegion()); tensorIt.GoToBegin(); int negatives = 0; while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iInitializeByItk( outputTensorImg.GetPointer() ); image->SetVolume( outputTensorImg->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dti_corrected"); SetDefaultNodeProperties(node, newname.toStdString()); nodes.push_back(node); // Corrected diffusion image // typedef itk::VectorImage ImageType; // ImageType::Pointer correctedVols = reconFilter->GetVectorImage(); // DiffusionImageType::Pointer correctedDiffusion = DiffusionImageType::New(); // correctedDiffusion->SetVectorImage(correctedVols); // correctedDiffusion->SetDirections(vols->GetDirections()); // correctedDiffusion->SetB_Value(vols->GetB_Value()); // correctedDiffusion->InitializeFromVectorImage(); // mitk::DataNode::Pointer diffNode = mitk::DataNode::New(); // diffNode->SetData( correctedDiffusion ); // QString diffname; // diffname = diffname.append(nodename.c_str()); // diffname = diffname.append("corrDiff"); // SetDefaultNodeProperties(diffNode, diffname.toStdString()); // nodes.push_back(diffNode); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); } } void QmitkTensorReconstructionView::ItkTensorReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // TENSOR RECONSTRUCTION clock.Start(); MITK_DEBUG << "Tensor reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::DiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientDirectionContainerType; GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = vols->GetDirections()->Begin(); it != vols->GetDirections()->End(); it++) { gradientContainerCopy->push_back(it.Value()); } tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, vols->GetVectorImage() ); tensorReconstructionFilter->SetBValue(vols->GetReferenceBValue()); tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreshold->value() ); tensorReconstructionFilter->Update(); clock.Stop(); MITK_DEBUG << "took " << clock.GetMean() << "s."; // TENSORS TO DATATREE mitk::TensorImage::Pointer image = mitk::TensorImage::New(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer tensorImage; tensorImage = tensorReconstructionFilter->GetOutput(); // Check the tensor for negative eigenvalues if(m_Controls->m_CheckNegativeEigenvalues->isChecked()) { typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetRequestedRegion()); tensorIt.GoToBegin(); while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; //typedef itk::Tensor TensorType2; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iSetDirection( vols->GetVectorImage()->GetDirection() ); image->InitializeByItk( tensorImage.GetPointer() ); image->SetVolume( tensorReconstructionFilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dti"); SetDefaultNodeProperties(node, newname.toStdString()); nodes.push_back(node); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return; } } void QmitkTensorReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name) { node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) ); node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) ); node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New()); node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New()); node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2)); node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1)); node->SetProperty( "visible", mitk::BoolProperty::New( true ) ); node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) ); node->SetProperty ("layer", mitk::IntProperty::New(100)); node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) ); node->SetProperty( "name", mitk::StringProperty::New(name) ); } void QmitkTensorReconstructionView::TensorsToDWI() { DoTensorsToDWI(m_TensorImages); } void QmitkTensorReconstructionView::TensorsToQbi() { for (unsigned int i=0; isize(); i++) { mitk::DataNode::Pointer tensorImageNode = m_TensorImages->at(i); MITK_INFO << "starting Q-Ball estimation"; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(dynamic_cast(tensorImageNode->GetData()), itkvol); typedef itk::TensorImageToQBallImageFilter< TTensorPixelType, TTensorPixelType > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->Update(); typedef itk::Vector OutputPixelType; typedef itk::Image OutputImageType; mitk::QBallImage::Pointer image = mitk::QBallImage::New(); OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(tensorImageNode->GetName().c_str()); newname = newname.append("_qbi"); node->SetName(newname.toAscii()); GetDefaultDataStorage()->Add(node); } } void QmitkTensorReconstructionView::OnSelectionChanged( std::vector nodes ) { m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); bool foundDwiVolume = false; bool foundTensorVolume = false; m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_DiffusionImage = NULL; m_TensorImage = NULL; m_Controls->m_InputData->setTitle("Please Select Input Data"); // iterate selection for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if (node.IsNull()) continue; // only look at interesting types if(dynamic_cast*>(node->GetData())) { foundDwiVolume = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_DiffusionImages->push_back(node); m_DiffusionImage = node; } else if(dynamic_cast(node->GetData())) { foundTensorVolume = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_TensorImages->push_back(node); m_TensorImage = node; } } m_Controls->m_StartReconstruction->setEnabled(foundDwiVolume); m_Controls->m_TensorsToDWIButton->setEnabled(foundTensorVolume); m_Controls->m_TensorsToQbiButton->setEnabled(foundTensorVolume); if (foundDwiVolume || foundTensorVolume) m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_ResidualButton->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PercentagesOfOutliers->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PerSliceView->setEnabled(foundDwiVolume && foundTensorVolume); } template -QmitkTensorReconstructionView::GradientListType::Pointer QmitkTensorReconstructionView::MakeGradientList() +itk::VectorContainer >::Pointer +QmitkTensorReconstructionView::MakeGradientList() { - QmitkTensorReconstructionView::GradientListType::Pointer retval = GradientListType::New(); + itk::VectorContainer >::Pointer retval = + itk::VectorContainer >::New(); vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); for(int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval->push_back(v); } // Add 0 vector for B0 - GradientType v(0.0); + vnl_vector_fixed v; + v.fill(0.0); retval->push_back(v); return retval; } void QmitkTensorReconstructionView::DoTensorsToDWI(mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { std::string nodename; (*itemiter)->GetStringProperty("name", nodename); mitk::TensorImage* vol = static_cast((*itemiter)->GetData()); ++itemiter; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(vol, itkvol); typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; - FilterType::GradientListType::Pointer gradientList; + FilterType::GradientListPointerType gradientList = FilterType::GradientListType::New(); switch(m_Controls->m_TensorsToDWINumDirsSelect->currentIndex()) { case 0: gradientList = MakeGradientList<12>(); break; case 1: gradientList = MakeGradientList<42>(); break; case 2: gradientList = MakeGradientList<92>(); break; case 3: gradientList = MakeGradientList<162>(); break; case 4: gradientList = MakeGradientList<252>(); break; case 5: gradientList = MakeGradientList<362>(); break; case 6: gradientList = MakeGradientList<492>(); break; case 7: gradientList = MakeGradientList<642>(); break; case 8: gradientList = MakeGradientList<812>(); break; case 9: gradientList = MakeGradientList<1002>(); break; default: gradientList = MakeGradientList<92>(); } double bVal = m_Controls->m_TensorsToDWIBValueEdit->text().toDouble(); // DWI ESTIMATION clock.Start(); MBI_INFO << "DWI Estimation "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "DWI Estimation for %s", nodename.c_str()).toAscii()); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->SetBValue(bVal); filter->SetGradientList(gradientList); //filter->SetNumberOfThreads(1); filter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMean() << "s."; // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetReferenceBValue(bVal); image->SetDirections(gradientList); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dwi"); node->SetName(newname.toAscii()); nodes.push_back(node); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "DWI estimation failed:", ex.GetDescription()); return ; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.h index 97671cea8a..cd927b796a 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.h @@ -1,115 +1,112 @@ /*=================================================================== 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 _QMITKTENSORRECONSTRUCTIONVIEW_H_INCLUDED #define _QMITKTENSORRECONSTRUCTIONVIEW_H_INCLUDED #include #include #include "ui_QmitkTensorReconstructionViewControls.h" #include #include typedef short DiffusionPixelType; struct TrSelListener; /*! * \ingroup org_mitk_gui_qt_tensorreconstruction_internal * * \brief QmitkTensorReconstructionView * * Document your class here. * * \sa QmitkFunctionality */ class QmitkTensorReconstructionView : public QmitkFunctionality { friend struct TrSelListener; // this is needed for all Qt objects that should have a MOC object (everything that derives from QObject) Q_OBJECT public: static const std::string VIEW_ID; QmitkTensorReconstructionView(); virtual ~QmitkTensorReconstructionView(); virtual void CreateQtPartControl(QWidget *parent); /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); /// \brief Called when the functionality is activated virtual void Activated(); virtual void Deactivated(); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); static const int nrconvkernels; protected slots: void TensorsToQbi(); void TensorsToDWI(); void DoTensorsToDWI(mitk::DataStorage::SetOfObjects::Pointer inImages); void Advanced1CheckboxClicked(); void Reconstruct(); void ResidualCalculation(); void ResidualClicked(int slice, int volume); protected: void ItkTensorReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages); void TeemTensorReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages); void TensorReconstructionWithCorr(mitk::DataStorage::SetOfObjects::Pointer inImages); void OnSelectionChanged( std::vector nodes ); Ui::QmitkTensorReconstructionViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; - typedef vnl_vector_fixed< double, 3 > GradientType; - typedef itk::VectorContainer< unsigned int, GradientType > GradientListType; - - template GradientListType::Pointer MakeGradientList(); + template itk::VectorContainer >::Pointer MakeGradientList(); template void TemplatedAnalyticalTensorReconstruction(mitk::DiffusionImage* vols, float lambda, std::string nodename, std::vector* nodes, int normalization); void SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name); mitk::DataNode::Pointer m_DiffusionImage; mitk::DataNode::Pointer m_TensorImage; mitk::DataStorage::SetOfObjects::Pointer m_DiffusionImages; mitk::DataStorage::SetOfObjects::Pointer m_TensorImages; }; #endif // _QMITKTENSORRECONSTRUCTIONVIEW_H_INCLUDED