diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp index c341de3251..db4cb6c8d7 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp @@ -1,314 +1,321 @@ /*=================================================================== 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 "mitkIOUtil.h" #include template< typename DiffusionPixelType> mitk::DWIHeadMotionCorrectionFilter - ::DWIHeadMotionCorrectionFilter(): -m_CurrentStep(0), +::DWIHeadMotionCorrectionFilter(): + m_CurrentStep(0), m_Steps(100), m_IsInValidState(true), - m_AbortRegistration(false) + 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 - registeredB0Image->SetImportVolume( registrationMethod->GetResampledMovingImage()->GetData(), - i, 0, mitk::Image::ReferenceMemory ); + // import volume to the inter-results + registeredB0Image->SetImportVolume( registrationMethod->GetResampledMovingImage()->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->GetB_ValueMap() ); + 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 registeredWeighted->SetImportVolume( b0referenceImage->GetData(), - 0,0, mitk::Image::CopyMemory ); + 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 registeredWeighted->SetImportVolume( weightedRegistrationMethod->GetResampledMovingImage()->GetData(), - i+1, 0, mitk::Image::CopyMemory); + 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->GetB_Value() ); 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()->SetB_Value(output->GetB_Value()); this->GetOutput()->SetDirections(output->GetDirections()); this->GetOutput()->SetMeasurementFrame(output->GetMeasurementFrame()); 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 69d9bc0ce5..272ef94c49 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h @@ -1,129 +1,131 @@ /*=================================================================== 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 ) itkNewMacro(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_Steps; unsigned long m_CurrentStep; 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) { // 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() ); } }; } //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 ec983dad04..c4ba37e886 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp @@ -1,140 +1,142 @@ /*=================================================================== 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_EstimatedParameters(NULL), + 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::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; } diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h index df66ef477e..256d29b007 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidImageRegistrationMethod.h @@ -1,452 +1,503 @@ /*=================================================================== 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 MITKPYRAMIDIMAGEREGISTRATION_H #define MITKPYRAMIDIMAGEREGISTRATION_H #include #include #include #include "mitkImage.h" #include "mitkBaseProcess.h" #include "mitkPyramidRegistrationMethodHelper.h" #include #include "mitkImageToItk.h" #include "mitkITKImageImport.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 DiffusionCore_EXPORT PyramidImageRegistrationMethod : public itk::Object { public: /** Typedefs */ mitkClassMacro(PyramidImageRegistrationMethod, itk::Object) /** Smart pointer support */ itkNewMacro(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; /** 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 ); 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; bool m_CrossModalityRegistration; bool m_UseAffineTransform; bool m_UseWindowedSincInterpolator; double* m_EstimatedParameters; + ParametersType m_InitialParameters; + /** Control the verbosity of the regsitistration output */ bool m_Verbose; 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::AffineTransform< double > AffineTransformType; typedef itk::Euler3DTransform< double > RigidTransformType; 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; } + 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 ); // 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); // 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 ); 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::WindowedSincInterpolateImageFunction< ImageType, 7> WindowedSincInterpolatorType; typename WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New(); typename mitk::ImageToItk< ImageType >::Pointer reference_image = mitk::ImageToItk< ImageType >::New(); reference_image->SetInput( this->m_FixedImage ); reference_image->Update(); 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 ); + //resampler->SetNumberOfThreads(1); + resampler->SetInput( itkImage ); resampler->SetTransform( base_transform ); resampler->SetReferenceImage( reference_image->GetOutput() ); 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 0fd024ad79..23d306c8c0 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkPyramidRegistrationMethodHelper.h @@ -1,122 +1,156 @@ /*=================================================================== 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))((float, dimension)), ((short, dimension))((float, 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; itkNewMacro( PyramidOptControlCommand ) 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: itkNewMacro( OptimizerIterationCommand ) 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 566c4806ca..75430dc8b4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkSplitDWImageFilter.h @@ -1,145 +1,152 @@ /*=================================================================== 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. */ itkNewMacro(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..0455285422 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.push_back(*listIt); + ++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 7143d2a49c..2e1b39d6ba 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h @@ -1,157 +1,157 @@ /*=================================================================== 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 TensorImageToDiffusionImageFilter Self; typedef ImageToImageFilter Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; itkTypeMacro (TensorImageToDiffusionImageFilter, ImageToImageFilter); itkStaticConstMacro (ImageDimension, unsigned int, OutputImageType::ImageDimension); itkNewMacro (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; } /** * @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; + 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..9670a723c8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx @@ -1,211 +1,217 @@ /*=================================================================== 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); 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); + unsigned int b0index = 0; + 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 ) + { + b0index = 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( 1.0 * b0 * exp ( -1.0 * bval * res ) ); } } } - out[0] = b0; + out[b0index] = 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/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp index a9b3452cae..ca87677a7d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp @@ -1,144 +1,146 @@ /*=================================================================== 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 directions and bvalue metaImage->SetDirections(this->m_GradientDirections); metaImage->SetB_Value(this->m_BValue); // set identity as measurement frame vnl_matrix_fixed< double, 3, 3 > measurement_frame; measurement_frame.set_identity(); metaImage->SetMeasurementFrame( measurement_frame ); 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()) ->SetB_Value(m_OutputCache->GetB_Value()); static_cast(this->GetOutput()) ->SetDirections(m_OutputCache->GetDirections()); static_cast(this->GetOutput()) ->SetMeasurementFrame(m_OutputCache->GetMeasurementFrame()); + 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 2b1323b4e3..2bd2f18497 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake @@ -1,9 +1,11 @@ set(MODULE_TESTS mitkFactoryRegistrationTest.cpp ) set(MODULE_CUSTOM_TESTS mitkPyramidImageRegistrationMethodTest.cpp mitkDWHeadMotionCorrectionTest.cpp + mitkConvertDWITypeTest.cpp + mitkExtractSingleShellTest.cpp ) diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp similarity index 81% copy from Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp copy to Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp index 5406afb055..fab2cf599b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp @@ -1,68 +1,63 @@ /*=================================================================== 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 "mitkDiffusionCoreObjectFactory.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."); RegisterDiffusionCoreObjectFactory(); 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( corrfilter->GetOutput() ); + dwiwriter->SetInput( dwimage ); 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/Testing/mitkDWHeadMotionCorrectionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp index 5406afb055..b81f80b328 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp @@ -1,68 +1,71 @@ /*=================================================================== 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 "mitkDiffusionCoreObjectFactory.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[] ) { MITK_TEST_BEGIN("mitkDWHeadMotionCorrectionTest"); MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output."); RegisterDiffusionCoreObjectFactory(); + 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(); mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::Pointer dwiwriter = mitk::NrrdDiffusionImageWriter< DiffusionPixelType >::New(); dwiwriter->SetInput( corrfilter->GetOutput() ); 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/Testing/mitkExtractSingleShellTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp new file mode 100644 index 0000000000..f6338063ce --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp @@ -0,0 +1,108 @@ +/*=================================================================== + +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 "mitkDiffusionCoreObjectFactory.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"); + + RegisterDiffusionCoreObjectFactory(); + + /* + 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; + int shellCounter = 0; + + 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 2264c3f07b..195b636b35 100644 --- a/Modules/DiffusionImaging/DiffusionCore/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake @@ -1,147 +1,149 @@ 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 IODataStructures/mitkDiffusionCoreObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp #IODataStructures/mitkRegistrationObject.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkPlanarFigureMapper3D.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/itkDwiGradientLengthCorrectionFilter.cpp + Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h # Registration Algorithms & Co. 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 Rendering/mitkPlanarFigureMapper3D.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/itkOdfMaximaExtractionFilter.h Algorithms/itkResampleDwiImageFilter.h Algorithms/itkDwiGradientLengthCorrectionFilter.h Algorithms/itkAdcImageFilter.h Algorithms/itkSplitDWImageFilter.h Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h Algorithms/mitkDiffusionImageToDiffusionImageFilter.h ) set( TOOL_FILES ) diff --git a/Modules/DiffusionImaging/FiberTracking/files.cmake b/Modules/DiffusionImaging/FiberTracking/files.cmake index 8ea747434e..b8546fdcd2 100644 --- a/Modules/DiffusionImaging/FiberTracking/files.cmake +++ b/Modules/DiffusionImaging/FiberTracking/files.cmake @@ -1,103 +1,104 @@ set(CPP_FILES # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkTrackvis.cpp # IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures IODataStructures/mitkFiberTrackingObjectFactory.cpp # Rendering Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp # Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp #Rendering/mitkPlanarFigureMapper3D.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 # Rendering Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h # Rendering/mitkFiberBundleXThreadMonitorMapper3D.h #Rendering/mitkPlanarFigureMapper3D.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h # IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h IODataStructures/FiberBundleX/mitkTrackvis.h IODataStructures/mitkFiberTrackingObjectFactory.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/mitkKspaceArtifact.h ) set(RESOURCE_FILES # Binary directory resources FiberTrackingLUTBaryCoords.bin FiberTrackingLUTIndices.bin # Shaders Shaders/mitkShaderFiberClipping.xml ) diff --git a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp index 30a0999123..081bfed7cb 100644 --- a/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp +++ b/Modules/ImageStatistics/mitkImageStatisticsCalculator.cpp @@ -1,1226 +1,1226 @@ /*=================================================================== 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 "mitkImageStatisticsCalculator.h" #include "mitkImageAccessByItk.h" #include "mitkImageCast.h" #include "mitkExtractImageFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #if ( ( VTK_MAJOR_VERSION <= 5 ) && ( VTK_MINOR_VERSION<=8) ) #include "mitkvtkLassoStencilSource.h" #else #include "vtkLassoStencilSource.h" #endif #include #include namespace mitk { ImageStatisticsCalculator::ImageStatisticsCalculator() : m_MaskingMode( MASKING_MODE_NONE ), m_MaskingModeChanged( false ), m_IgnorePixelValue(0.0), m_DoIgnorePixelValue(false), m_IgnorePixelValueChanged(false), m_PlanarFigureAxis (0), m_PlanarFigureSlice (0), m_PlanarFigureCoordinate0 (0), m_PlanarFigureCoordinate1 (0) { m_EmptyHistogram = HistogramType::New(); m_EmptyHistogram->SetMeasurementVectorSize(1); HistogramType::SizeType histogramSize(1); histogramSize.Fill( 256 ); m_EmptyHistogram->Initialize( histogramSize ); m_EmptyStatistics.Reset(); } ImageStatisticsCalculator::~ImageStatisticsCalculator() { } void ImageStatisticsCalculator::SetImage( const mitk::Image *image ) { if ( m_Image != image ) { m_Image = image; this->Modified(); unsigned int numberOfTimeSteps = image->GetTimeSteps(); // Initialize vectors to time-size of this image m_ImageHistogramVector.resize( numberOfTimeSteps ); m_MaskedImageHistogramVector.resize( numberOfTimeSteps ); m_PlanarFigureHistogramVector.resize( numberOfTimeSteps ); m_ImageStatisticsVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsVector.resize( numberOfTimeSteps ); m_ImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsTimeStampVector.resize( numberOfTimeSteps ); m_ImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_MaskedImageStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); m_PlanarFigureStatisticsCalculationTriggerVector.resize( numberOfTimeSteps ); for ( unsigned int t = 0; t < image->GetTimeSteps(); ++t ) { m_ImageStatisticsTimeStampVector[t].Modified(); m_ImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetImageMask( const mitk::Image *imageMask ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) { itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); } if ( m_ImageMask != imageMask ) { m_ImageMask = imageMask; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_MaskedImageStatisticsTimeStampVector[t].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetPlanarFigure( mitk::PlanarFigure *planarFigure ) { if ( m_Image.IsNull() ) { itkExceptionMacro( << "Image needs to be set first!" ); } if ( m_PlanarFigure != planarFigure ) { m_PlanarFigure = planarFigure; this->Modified(); for ( unsigned int t = 0; t < m_Image->GetTimeSteps(); ++t ) { m_PlanarFigureStatisticsTimeStampVector[t].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[t] = true; } } } void ImageStatisticsCalculator::SetMaskingMode( unsigned int mode ) { if ( m_MaskingMode != mode ) { m_MaskingMode = mode; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToNone() { if ( m_MaskingMode != MASKING_MODE_NONE ) { m_MaskingMode = MASKING_MODE_NONE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToImage() { if ( m_MaskingMode != MASKING_MODE_IMAGE ) { m_MaskingMode = MASKING_MODE_IMAGE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetMaskingModeToPlanarFigure() { if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) { m_MaskingMode = MASKING_MODE_PLANARFIGURE; m_MaskingModeChanged = true; this->Modified(); } } void ImageStatisticsCalculator::SetIgnorePixelValue(double value) { if ( m_IgnorePixelValue != value ) { m_IgnorePixelValue = value; if(m_DoIgnorePixelValue) { m_IgnorePixelValueChanged = true; } this->Modified(); } } double ImageStatisticsCalculator::GetIgnorePixelValue() { return m_IgnorePixelValue; } void ImageStatisticsCalculator::SetDoIgnorePixelValue(bool value) { if ( m_DoIgnorePixelValue != value ) { m_DoIgnorePixelValue = value; m_IgnorePixelValueChanged = true; this->Modified(); } } bool ImageStatisticsCalculator::GetDoIgnorePixelValue() { return m_DoIgnorePixelValue; } bool ImageStatisticsCalculator::ComputeStatistics( unsigned int timeStep ) { if (m_Image.IsNull() ) { mitkThrow() << "Image not set!"; } if (!m_Image->IsInitialized()) { mitkThrow() << "Image not initialized!"; } if ( m_Image->GetReferenceCount() == 1 ) { // Image no longer valid; we are the only ones to still hold a reference on it return false; } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } // If a mask was set but we are the only ones to still hold a reference on // it, delete it. if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) { m_ImageMask = NULL; } // Check if statistics is already up-to-date unsigned long imageMTime = m_ImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStampVector[timeStep].GetMTime(); unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStampVector[timeStep].GetMTime(); bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerVector[timeStep]; bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerVector[timeStep]; bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerVector[timeStep]; if ( !m_IgnorePixelValueChanged && ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_IMAGE) || (maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) { // Statistics is up to date! if ( m_MaskingModeChanged ) { m_MaskingModeChanged = false; return true; } else { return false; } } // Reset state changed flag m_MaskingModeChanged = false; m_IgnorePixelValueChanged = false; // Depending on masking mode, extract and/or generate the required image // and mask data from the user input this->ExtractImageAndMask( timeStep ); StatisticsContainer *statisticsContainer; HistogramContainer *histogramContainer; switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: if(!m_DoIgnorePixelValue) { statisticsContainer = &m_ImageStatisticsVector[timeStep]; histogramContainer = &m_ImageHistogramVector[timeStep]; m_ImageStatisticsTimeStampVector[timeStep].Modified(); m_ImageStatisticsCalculationTriggerVector[timeStep] = false; } else { statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; } break; case MASKING_MODE_IMAGE: statisticsContainer = &m_MaskedImageStatisticsVector[timeStep]; histogramContainer = &m_MaskedImageHistogramVector[timeStep]; m_MaskedImageStatisticsTimeStampVector[timeStep].Modified(); m_MaskedImageStatisticsCalculationTriggerVector[timeStep] = false; break; case MASKING_MODE_PLANARFIGURE: statisticsContainer = &m_PlanarFigureStatisticsVector[timeStep]; histogramContainer = &m_PlanarFigureHistogramVector[timeStep]; m_PlanarFigureStatisticsTimeStampVector[timeStep].Modified(); m_PlanarFigureStatisticsCalculationTriggerVector[timeStep] = false; break; } // Calculate statistics and histogram(s) if ( m_InternalImage->GetDimension() == 3 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 3, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 3, m_InternalImageMask3D.GetPointer(), statisticsContainer, histogramContainer ); } } else if ( m_InternalImage->GetDimension() == 2 ) { if ( m_MaskingMode == MASKING_MODE_NONE && !m_DoIgnorePixelValue ) { AccessFixedDimensionByItk_2( m_InternalImage, InternalCalculateStatisticsUnmasked, 2, statisticsContainer, histogramContainer ); } else { AccessFixedDimensionByItk_3( m_InternalImage, InternalCalculateStatisticsMasked, 2, m_InternalImageMask2D.GetPointer(), statisticsContainer, histogramContainer ); } } else { MITK_ERROR << "ImageStatistics: Image dimension not supported!"; } // Release unused image smart pointers to free memory m_InternalImage = mitk::Image::ConstPointer(); m_InternalImageMask3D = MaskImage3DType::Pointer(); m_InternalImageMask2D = MaskImage2DType::Pointer(); return true; } const ImageStatisticsCalculator::HistogramType * ImageStatisticsCalculator::GetHistogram( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return NULL; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep][label]; return m_ImageHistogramVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep][label]; } } const ImageStatisticsCalculator::HistogramContainer & ImageStatisticsCalculator::GetHistogramVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyHistogramContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageHistogramVector[timeStep]; return m_ImageHistogramVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageHistogramVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureHistogramVector[timeStep]; } } const ImageStatisticsCalculator::Statistics & ImageStatisticsCalculator::GetStatistics( unsigned int timeStep, unsigned int label ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatistics; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep][label]; return m_ImageStatisticsVector[timeStep][label]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep][label]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep][label]; } } const ImageStatisticsCalculator::StatisticsContainer & ImageStatisticsCalculator::GetStatisticsVector( unsigned int timeStep ) const { if ( m_Image.IsNull() || (timeStep >= m_Image->GetTimeSteps()) ) { return m_EmptyStatisticsContainer; } switch ( m_MaskingMode ) { case MASKING_MODE_NONE: default: { if(m_DoIgnorePixelValue) return m_MaskedImageStatisticsVector[timeStep]; return m_ImageStatisticsVector[timeStep]; } case MASKING_MODE_IMAGE: return m_MaskedImageStatisticsVector[timeStep]; case MASKING_MODE_PLANARFIGURE: return m_PlanarFigureStatisticsVector[timeStep]; } } void ImageStatisticsCalculator::ExtractImageAndMask( unsigned int timeStep ) { if ( m_Image.IsNull() ) { throw std::runtime_error( "Error: image empty!" ); } if ( timeStep >= m_Image->GetTimeSteps() ) { throw std::runtime_error( "Error: invalid time step!" ); } ImageTimeSelector::Pointer imageTimeSelector = ImageTimeSelector::New(); imageTimeSelector->SetInput( m_Image ); imageTimeSelector->SetTimeNr( timeStep ); imageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceImage = imageTimeSelector->GetOutput(); switch ( m_MaskingMode ) { case MASKING_MODE_NONE: { m_InternalImage = timeSliceImage; m_InternalImageMask2D = NULL; m_InternalImageMask3D = NULL; if(m_DoIgnorePixelValue) { if( m_InternalImage->GetDimension() == 3 ) { CastToItkImage( timeSliceImage, m_InternalImageMask3D ); m_InternalImageMask3D->FillBuffer(1); } if( m_InternalImage->GetDimension() == 2 ) { CastToItkImage( timeSliceImage, m_InternalImageMask2D ); m_InternalImageMask2D->FillBuffer(1); } } break; } case MASKING_MODE_IMAGE: { if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) { if ( timeStep < m_ImageMask->GetTimeSteps() ) { ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); maskedImageTimeSelector->SetInput( m_ImageMask ); maskedImageTimeSelector->SetTimeNr( timeStep ); maskedImageTimeSelector->UpdateLargestPossibleRegion(); mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); m_InternalImage = timeSliceImage; CastToItkImage( timeSliceMaskedImage, m_InternalImageMask3D ); } else { throw std::runtime_error( "Error: image mask has not enough time steps!" ); } } else { throw std::runtime_error( "Error: image mask empty!" ); } break; } case MASKING_MODE_PLANARFIGURE: { m_InternalImageMask2D = NULL; if ( m_PlanarFigure.IsNull() ) { throw std::runtime_error( "Error: planar figure empty!" ); } if ( !m_PlanarFigure->IsClosed() ) { throw std::runtime_error( "Masking not possible for non-closed figures" ); } const Geometry3D *imageGeometry = timeSliceImage->GetGeometry(); if ( imageGeometry == NULL ) { throw std::runtime_error( "Image geometry invalid!" ); } const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); if ( planarFigureGeometry2D == NULL ) { throw std::runtime_error( "Planar-Figure not yet initialized!" ); } const PlaneGeometry *planarFigureGeometry = dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); if ( planarFigureGeometry == NULL ) { throw std::runtime_error( "Non-planar planar figures not supported!" ); } // Find principal direction of PlanarFigure in input image unsigned int axis; if ( !this->GetPrincipalAxis( imageGeometry, planarFigureGeometry->GetNormal(), axis ) ) { throw std::runtime_error( "Non-aligned planar figures not supported!" ); } m_PlanarFigureAxis = axis; // Find slice number corresponding to PlanarFigure in input image MaskImage3DType::IndexType index; imageGeometry->WorldToIndex( planarFigureGeometry->GetOrigin(), index ); unsigned int slice = index[axis]; m_PlanarFigureSlice = slice; // Extract slice with given position and direction from image unsigned int dimension = timeSliceImage->GetDimension(); if (dimension != 2) { ExtractImageFilter::Pointer imageExtractor = ExtractImageFilter::New(); imageExtractor->SetInput( timeSliceImage ); imageExtractor->SetSliceDimension( axis ); imageExtractor->SetSliceIndex( slice ); imageExtractor->Update(); m_InternalImage = imageExtractor->GetOutput(); } else { m_InternalImage = timeSliceImage; } // Compute mask from PlanarFigure AccessFixedDimensionByItk_1( m_InternalImage, InternalCalculateMaskFromPlanarFigure, 2, axis ); } } if(m_DoIgnorePixelValue) { if ( m_InternalImage->GetDimension() == 3 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 3, m_InternalImageMask3D.GetPointer() ); } else if ( m_InternalImage->GetDimension() == 2 ) { AccessFixedDimensionByItk_1( m_InternalImage, InternalMaskIgnoredPixels, 2, m_InternalImageMask2D.GetPointer() ); } } } bool ImageStatisticsCalculator::GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, unsigned int &axis ) { vector.Normalize(); for ( unsigned int i = 0; i < 3; ++i ) { Vector3D axisVector = geometry->GetAxisVector( i ); axisVector.Normalize(); if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) { axis = i; return true; } } return false; } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsUnmasked( const itk::Image< TPixel, VImageDimension > *image, StatisticsContainer *statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > HistogramGeneratorType; statisticsContainer->clear(); histogramContainer->clear(); // Progress listening... typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate ); // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator // does not (yet?) support progress reporting this->InvokeEvent( itk::StartEvent() ); for ( unsigned int i = 0; i < 100; ++i ) { this->UnmaskedStatisticsProgressUpdate(); } // Calculate statistics (separate filter) typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( image ); unsigned long observerTag = statisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); statisticsFilter->Update(); statisticsFilter->RemoveObserver( observerTag ); this->InvokeEvent( itk::EndEvent() ); // Calculate minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( image ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); Statistics statistics; statistics.Reset(); statistics.Label = 1; statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); statistics.Min = statisticsFilter->GetMinimum(); statistics.Max = statisticsFilter->GetMaximum(); statistics.Mean = statisticsFilter->GetMean(); statistics.Median = 0.0; statistics.Sigma = statisticsFilter->GetSigma(); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); statistics.MinIndex.set_size(image->GetImageDimension()); statistics.MaxIndex.set_size(image->GetImageDimension()); - for (int i=0; iGetIndexOfMaximum()[i]; statistics.MinIndex[i] = minMaxFilter->GetIndexOfMinimum()[i]; } statisticsContainer->push_back( statistics ); // Calculate histogram typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); histogramGenerator->SetInput( image ); histogramGenerator->SetMarginalScale( 100 ); histogramGenerator->SetNumberOfBins( 768 ); histogramGenerator->SetHistogramMin( statistics.Min ); histogramGenerator->SetHistogramMax( statistics.Max ); histogramGenerator->Compute(); histogramContainer->push_back( histogramGenerator->GetOutput() ); } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalMaskIgnoredPixels( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; itk::ImageRegionIterator itmask(maskImage, maskImage->GetLargestPossibleRegion()); itk::ImageRegionConstIterator itimage(image, image->GetLargestPossibleRegion()); itmask.GoToBegin(); itimage.GoToBegin(); while( !itmask.IsAtEnd() ) { if(m_IgnorePixelValue == itimage.Get()) { itmask.Set(0); } ++itmask; ++itimage; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateStatisticsMasked( const itk::Image< TPixel, VImageDimension > *image, itk::Image< unsigned short, VImageDimension > *maskImage, StatisticsContainer* statisticsContainer, HistogramContainer* histogramContainer ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::Image< unsigned short, VImageDimension > MaskImageType; typedef typename ImageType::IndexType IndexType; typedef typename ImageType::PointType PointType; typedef typename ImageType::SpacingType SpacingType; typedef itk::LabelStatisticsImageFilter< ImageType, MaskImageType > LabelStatisticsFilterType; typedef itk::ChangeInformationImageFilter< MaskImageType > ChangeInformationFilterType; typedef itk::ExtractImageFilter< ImageType, ImageType > ExtractImageFilterType; statisticsContainer->clear(); histogramContainer->clear(); // Make sure that mask is set if ( maskImage == NULL ) { itkExceptionMacro( << "Mask image needs to be set!" ); } // Make sure that spacing of mask and image are the same SpacingType imageSpacing = image->GetSpacing(); SpacingType maskSpacing = maskImage->GetSpacing(); PointType zeroPoint; zeroPoint.Fill( 0.0 ); if ( (zeroPoint + imageSpacing).SquaredEuclideanDistanceTo( (zeroPoint + maskSpacing) ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same spacing as image! (Image spacing: " << imageSpacing << "; Mask spacing: " << maskSpacing << ")" ); } // Make sure that orientation of mask and image are the same typedef typename ImageType::DirectionType DirectionType; DirectionType imageDirection = image->GetDirection(); DirectionType maskDirection = maskImage->GetDirection(); for( int i = 0; i < imageDirection.ColumnDimensions; ++i ) { for( int j = 0; j < imageDirection.ColumnDimensions; ++j ) { double differenceDirection = imageDirection[i][j] - maskDirection[i][j]; if ( fabs( differenceDirection ) > mitk::eps ) { itkExceptionMacro( << "Mask needs to have same direction as image! (Image direction: " << imageDirection << "; Mask direction: " << maskDirection << ")" ); } } } // Make sure that the voxels of mask and image are correctly "aligned", i.e., voxel boundaries are the same in both images PointType imageOrigin = image->GetOrigin(); PointType maskOrigin = maskImage->GetOrigin(); long offset[ImageType::ImageDimension]; typedef itk::ContinuousIndex ContinousIndexType; ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; image->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); image->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); if ( fabs( misalignment ) > mitk::eps ) { itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); } double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; offset[i] = (int) indexCoordDistance + image->GetBufferedRegion().GetIndex()[i]; } // Adapt the origin and region (index/size) of the mask so that the origin of both are the same typename ChangeInformationFilterType::Pointer adaptMaskFilter; adaptMaskFilter = ChangeInformationFilterType::New(); adaptMaskFilter->ChangeOriginOn(); adaptMaskFilter->ChangeRegionOn(); adaptMaskFilter->SetInput( maskImage ); adaptMaskFilter->SetOutputOrigin( image->GetOrigin() ); adaptMaskFilter->SetOutputOffset( offset ); adaptMaskFilter->Update(); typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); // Make sure that mask region is contained within image region if ( !image->GetLargestPossibleRegion().IsInside( adaptedMaskImage->GetLargestPossibleRegion() ) ) { itkExceptionMacro( << "Mask region needs to be inside of image region! (Image region: " << image->GetLargestPossibleRegion() << "; Mask region: " << adaptedMaskImage->GetLargestPossibleRegion() << ")" ); } // If mask region is smaller than image region, extract the sub-sampled region from the original image typename ImageType::SizeType imageSize = image->GetBufferedRegion().GetSize(); typename ImageType::SizeType maskSize = maskImage->GetBufferedRegion().GetSize(); bool maskSmallerImage = false; for ( unsigned int i = 0; i < ImageType::ImageDimension; ++i ) { if ( maskSize[i] < imageSize[i] ) { maskSmallerImage = true; } } typename ImageType::ConstPointer adaptedImage; if ( maskSmallerImage ) { typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); extractImageFilter->SetInput( image ); extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); extractImageFilter->Update(); adaptedImage = extractImageFilter->GetOutput(); } else { adaptedImage = image; } // Initialize Filter typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); statisticsFilter->SetInput( adaptedImage ); statisticsFilter->Update(); int numberOfBins = ( m_DoIgnorePixelValue && (m_MaskingMode == MASKING_MODE_NONE) ) ? 768 : 384; typename LabelStatisticsFilterType::Pointer labelStatisticsFilter; labelStatisticsFilter = LabelStatisticsFilterType::New(); labelStatisticsFilter->SetInput( adaptedImage ); labelStatisticsFilter->SetLabelInput( adaptedMaskImage ); labelStatisticsFilter->UseHistogramsOn(); labelStatisticsFilter->SetHistogramParameters( numberOfBins, statisticsFilter->GetMinimum(), statisticsFilter->GetMaximum() ); // Add progress listening typedef itk::SimpleMemberCommand< ImageStatisticsCalculator > ITKCommandType; ITKCommandType::Pointer progressListener; progressListener = ITKCommandType::New(); progressListener->SetCallbackFunction( this, &ImageStatisticsCalculator::MaskedStatisticsProgressUpdate ); unsigned long observerTag = labelStatisticsFilter->AddObserver( itk::ProgressEvent(), progressListener ); // Execute filter this->InvokeEvent( itk::StartEvent() ); // Make sure that only the mask region is considered (otherwise, if the mask region is smaller // than the image region, the Update() would result in an exception). labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); // Execute the filter labelStatisticsFilter->Update(); this->InvokeEvent( itk::EndEvent() ); labelStatisticsFilter->RemoveObserver( observerTag ); // Find all relevant labels of mask (other than 0) std::list< int > relevantLabels; bool maskNonEmpty = false; unsigned int i; for ( i = 1; i < 4096; ++i ) { if ( labelStatisticsFilter->HasLabel( i ) ) { relevantLabels.push_back( i ); maskNonEmpty = true; } } if ( maskNonEmpty ) { std::list< int >::iterator it; for ( it = relevantLabels.begin(), i = 0; it != relevantLabels.end(); ++it, ++i ) { histogramContainer->push_back( HistogramType::ConstPointer( labelStatisticsFilter->GetHistogram( (*it) ) ) ); Statistics statistics; statistics.Label = (*it); statistics.N = labelStatisticsFilter->GetCount( *it ); statistics.Min = labelStatisticsFilter->GetMinimum( *it ); statistics.Max = labelStatisticsFilter->GetMaximum( *it ); statistics.Mean = labelStatisticsFilter->GetMean( *it ); statistics.Median = labelStatisticsFilter->GetMedian( *it ); statistics.Sigma = labelStatisticsFilter->GetSigma( *it ); statistics.RMS = sqrt( statistics.Mean * statistics.Mean + statistics.Sigma * statistics.Sigma ); // restrict image to mask area for min/max index calculation typedef itk::MaskImageFilter< ImageType, MaskImageType, ImageType > MaskImageFilterType; typename MaskImageFilterType::Pointer masker = MaskImageFilterType::New(); masker->SetOutsideValue( (statistics.Min+statistics.Max)/2 ); masker->SetInput1(adaptedImage); masker->SetInput2(adaptedMaskImage); masker->Update(); // get index of minimum and maximum typedef itk::MinimumMaximumImageCalculator< ImageType > MinMaxFilterType; typename MinMaxFilterType::Pointer minMaxFilter = MinMaxFilterType::New(); minMaxFilter->SetImage( masker->GetOutput() ); unsigned long observerTag2 = minMaxFilter->AddObserver( itk::ProgressEvent(), progressListener ); minMaxFilter->Compute(); minMaxFilter->RemoveObserver( observerTag2 ); this->InvokeEvent( itk::EndEvent() ); statistics.MinIndex.set_size(adaptedImage->GetImageDimension()); statistics.MaxIndex.set_size(adaptedImage->GetImageDimension()); typename MinMaxFilterType::IndexType tempMaxIndex = minMaxFilter->GetIndexOfMaximum(); typename MinMaxFilterType::IndexType tempMinIndex = minMaxFilter->GetIndexOfMinimum(); // FIX BUG 14644 //If a PlanarFigure is used for segmentation the //adaptedImage is a single slice (2D). Adding the // 3. dimension. if (m_MaskingMode == MASKING_MODE_PLANARFIGURE && m_Image->GetDimension()==3) { statistics.MaxIndex.set_size(m_Image->GetDimension()); statistics.MaxIndex[m_PlanarFigureCoordinate0]=tempMaxIndex[0]; statistics.MaxIndex[m_PlanarFigureCoordinate1]=tempMaxIndex[1]; statistics.MaxIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; statistics.MinIndex.set_size(m_Image->GetDimension()); statistics.MinIndex[m_PlanarFigureCoordinate0]=tempMinIndex[0]; statistics.MinIndex[m_PlanarFigureCoordinate1]=tempMinIndex[1]; statistics.MinIndex[m_PlanarFigureAxis]=m_PlanarFigureSlice; } else { for (int i = 0; ipush_back( statistics ); } } else { histogramContainer->push_back( HistogramType::ConstPointer( m_EmptyHistogram ) ); statisticsContainer->push_back( Statistics() );; } } template < typename TPixel, unsigned int VImageDimension > void ImageStatisticsCalculator::InternalCalculateMaskFromPlanarFigure( const itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) { typedef itk::Image< TPixel, VImageDimension > ImageType; typedef itk::CastImageFilter< ImageType, MaskImage2DType > CastFilterType; // Generate mask image as new image with same header as input image and // initialize with "1". typename CastFilterType::Pointer castFilter = CastFilterType::New(); castFilter->SetInput( image ); castFilter->Update(); castFilter->GetOutput()->FillBuffer( 1 ); // all PolylinePoints of the PlanarFigure are stored in a vtkPoints object. // These points are used by the vtkLassoStencilSource to create // a vtkImageStencil. const mitk::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); const mitk::Geometry3D *imageGeometry3D = m_Image->GetGeometry( 0 ); // Determine x- and y-dimensions depending on principal axis int i0, i1; switch ( axis ) { case 0: i0 = 1; i1 = 2; break; case 1: i0 = 0; i1 = 2; break; case 2: default: i0 = 0; i1 = 1; break; } m_PlanarFigureCoordinate0= i0; m_PlanarFigureCoordinate1= i1; // store the polyline contour as vtkPoints object bool outOfBounds = false; vtkSmartPointer points = vtkSmartPointer::New(); typename PlanarFigure::PolyLineType::const_iterator it; for ( it = planarFigurePolyline.begin(); it != planarFigurePolyline.end(); ++it ) { Point3D point3D; // Convert 2D point back to the local index coordinates of the selected // image planarFigureGeometry2D->Map( it->Point, point3D ); // Polygons (partially) outside of the image bounds can not be processed // further due to a bug in vtkPolyDataToImageStencil if ( !imageGeometry3D->IsInside( point3D ) ) { outOfBounds = true; } imageGeometry3D->WorldToIndex( point3D, point3D ); points->InsertNextPoint( point3D[i0], point3D[i1], 0 ); } // mark a malformed 2D planar figure ( i.e. area = 0 ) as out of bounds // this can happen when all control points of a rectangle lie on the same line = two of the three extents are zero double bounds[6] = {0, 0, 0, 0, 0, 0}; points->GetBounds( bounds ); bool extent_x = (fabs(bounds[0] - bounds[1])) < mitk::eps; bool extent_y = (fabs(bounds[2] - bounds[3])) < mitk::eps; bool extent_z = (fabs(bounds[4] - bounds[5])) < mitk::eps; // throw an exception if a closed planar figure is deformed, i.e. has only one non-zero extent if ( m_PlanarFigure->IsClosed() && ((extent_x && extent_y) || (extent_x && extent_z) || (extent_y && extent_z))) { mitkThrow() << "Figure has a zero area and cannot be used for masking."; } if ( outOfBounds ) { throw std::runtime_error( "Figure at least partially outside of image bounds!" ); } // create a vtkLassoStencilSource and set the points of the Polygon vtkSmartPointer lassoStencil = vtkSmartPointer::New(); lassoStencil->SetShapeToPolygon(); lassoStencil->SetPoints( points ); // Export from ITK to VTK (to use a VTK filter) typedef itk::VTKImageImport< MaskImage2DType > ImageImportType; typedef itk::VTKImageExport< MaskImage2DType > ImageExportType; typename ImageExportType::Pointer itkExporter = ImageExportType::New(); itkExporter->SetInput( castFilter->GetOutput() ); vtkSmartPointer vtkImporter = vtkSmartPointer::New(); this->ConnectPipelines( itkExporter, vtkImporter ); // Apply the generated image stencil to the input image vtkSmartPointer imageStencilFilter = vtkSmartPointer::New(); imageStencilFilter->SetInputConnection( vtkImporter->GetOutputPort() ); imageStencilFilter->SetStencil( lassoStencil->GetOutput() ); imageStencilFilter->ReverseStencilOff(); imageStencilFilter->SetBackgroundValue( 0 ); imageStencilFilter->Update(); // Export from VTK back to ITK vtkSmartPointer vtkExporter = vtkImageExport::New(); // TODO: this is WRONG, should be vtkSmartPointer::New(), but bug # 14455 //vtkSmartPointer vtkExporter = vtkSmartPointer::New(); vtkExporter->SetInputConnection( imageStencilFilter->GetOutputPort() ); vtkExporter->Update(); typename ImageImportType::Pointer itkImporter = ImageImportType::New(); this->ConnectPipelines( vtkExporter, itkImporter ); itkImporter->Update(); // Store mask m_InternalImageMask2D = itkImporter->GetOutput(); } void ImageStatisticsCalculator::UnmaskedStatisticsProgressUpdate() { // Need to throw away every second progress event to reach a final count of // 100 since two consecutive filters are used in this case static int updateCounter = 0; if ( updateCounter++ % 2 == 0 ) { this->InvokeEvent( itk::ProgressEvent() ); } } void ImageStatisticsCalculator::MaskedStatisticsProgressUpdate() { this->InvokeEvent( itk::ProgressEvent() ); } } 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 ae26054dc0..df6f1c93df 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp @@ -1,1031 +1,1034 @@ /*=================================================================== 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(const QmitkTensorReconstructionView& other) { Q_UNUSED(other) throw std::runtime_error("Copy constructor not implemented"); } 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(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->GetB_Value()); - filter->SetGradientList(gradients); + filter->SetGradientList( diffImage->GetDirections() ); 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->SetB_Value(diffImage->GetB_Value()); - image->SetDirections(gradients); + image->SetDirections(diffImage->GetDirections()); 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); std::vector b0Indices = image->GetB0Indices(); 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(); int size = lookupTable->GetTable()->GetSize(); 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->GetB_Value()); 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->GetB_Value()); tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreshold->value() ); tensorReconstructionFilter->Update(); clock.Stop(); MITK_DEBUG << "took " << clock.GetMeanTime() << "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 (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.GetMeanTime() << "s."; // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(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 d74669300e..c86673eea5 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.h @@ -1,116 +1,113 @@ /*=================================================================== 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(); QmitkTensorReconstructionView(const QmitkTensorReconstructionView& other); 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