diff --git a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h index 9c9357be60..3533d75b12 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h @@ -1,96 +1,98 @@ /*=================================================================== 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 MITKDIFFUSIONIMAGECORRECTIONFILTER_H #define MITKDIFFUSIONIMAGECORRECTIONFILTER_H #include "mitkImageSource.h" #include namespace mitk { /** * @class DiffusionImageCorrectionFilter */ class MITKDIFFUSIONCORE_EXPORT DiffusionImageCorrectionFilter : public ImageSource { public: /** class macros */ mitkClassMacro( DiffusionImageCorrectionFilter, ImageSource ) itkSimpleNewMacro(Self) typedef short DiffusionPixelType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef vnl_matrix_fixed< double, 3, 3 > TransformMatrixType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef GradientDirectionContainerType::Pointer GradientDirectionContainerPointerType; typedef std::vector< TransformMatrixType > TransformsVectorType; typedef Superclass::OutputType DiffusionImageType; typedef DiffusionImageType::Pointer DiffusionImageTypePointer; typedef itk::VectorImage ImageType; /** * @brief Set the mitk image ( a 3d+t image ) which is to be reinterpreted as dw image * @param mitkImage */ void SetImage( DiffusionImageTypePointer input ) { m_SourceImage = input; } /** * @brief Correct each gradient direction according to the given transform * * The size of the input is expected to correspond to the count of gradient images in the image. */ void CorrectDirections( const TransformsVectorType& ); /** * @brief Correct all gradient directions according to the given transform * * This will apply the same rotation to all directions. */ void CorrectDirections( const TransformMatrixType& ); + virtual void GenerateOutputInformation() override {} + protected: DiffusionImageCorrectionFilter(); virtual ~DiffusionImageCorrectionFilter() {} /** * @brief Get the rotation component following the Finite Strain * * For a given transformation \f$A\f$ its rotation component is defined as \f$ (AA^{T})^{-1/2}\f$. * * The computation first computes \f$ B = AA^T \f$ and then estimates the square root. Square root of * diagonal matrices is defined as * \f$ S = Q * \sqrt{C} * Q^{-1} \f$ with \f$ C \f$ having the eigenvalues on the diagonal. * */ TransformMatrixType GetRotationComponent( const TransformMatrixType& ); DiffusionImageTypePointer m_SourceImage; }; } #endif // MITKDIFFUSIONIMAGECORRECTIONFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.h index 1eaac221d6..0a104c7fe2 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.h @@ -1,107 +1,111 @@ /*=================================================================== 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 DIFFUSIONIMAGECREATIONFILTER_H #define DIFFUSIONIMAGECREATIONFILTER_H #include "mitkImageToImageFilter.h" #include "mitkDiffusionPropertyHelper.h" #include namespace mitk { /** * @brief The DiffusionImageHeaderDescriptor struct bundles the necessary diffusion-weighted image header meta information. */ struct MITKDIFFUSIONCORE_EXPORT DiffusionImageHeaderDescriptor { DiffusionImageHeaderDescriptor() : m_GradientDirections(nullptr) { m_BValue = -1; } float m_BValue; mitk::DiffusionPropertyHelper::BValueMapType m_BValueMapType; mitk::DiffusionPropertyHelper::MeasurementFrameType m_MeasurementFrame; mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer m_GradientDirections; }; /** * @brief The DiffusionImageCreationFilter class creates a diffusion-weighted image out of a * given 3D+t regular image and sufficient additional information about gradient directions and b values * * For easier use, the filter supports the usage of reference images. Here a diffusion-weighted (dw) image is expected and the meta * information of this image will be used for the output dw image. The diffusion information can also be specified by setting the HeaderDescriptor * directly (SetHeaderDescriptor). * * At least one of reference image or descriptor must be used, otherwise an exception is thrown. */ class MITKDIFFUSIONCORE_EXPORT DiffusionImageCreationFilter : public ImageToImageFilter { public: /** class macros */ mitkClassMacro( DiffusionImageCreationFilter, ImageToImageFilter ) itkNewMacro(Self) typedef short DiffusionPixelType; typedef mitk::DiffusionPropertyHelper DPH; typedef mitk::Image OutputType; typedef mitk::DiffusionPropertyHelper::ImageType VectorImageType; typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; typedef mitk::DiffusionPropertyHelper::MeasurementFrameType MeasurementFrameType; typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; /** * @brief SetReferenceImage Set a diffusion image as reference, i.e. the header information will be extracted from it * @param reference_image A reference diffusion-weigted image - will throw exception of the input is not DW image */ void SetReferenceImage( mitk::Image::Pointer reference_image ); /** * @brief SetHeaderDescriptor set the information to be used with the dw image * @param header_descriptor * * \sa DiffusionImageHeaderDescriptor */ void SetHeaderDescriptor( DiffusionImageHeaderDescriptor header_descriptor ); virtual void GenerateData() override; + virtual void GenerateOutputInformation() override; + protected: DiffusionImageCreationFilter(); ~DiffusionImageCreationFilter(); GradientDirectionContainerType::Pointer InternalGetGradientDirections( ); MeasurementFrameType InternalGetMeasurementFrame(); float InternalGetBValue(); - VectorImageType::Pointer RemapIntoVectorImage( mitk::Image::Pointer input); + mitk::Image::Pointer RemapIntoVectorImage( mitk::Image::Pointer input); mitk::Image::Pointer m_ReferenceImage; + mitk::Image::Pointer m_OutputCache; + DiffusionImageHeaderDescriptor m_HeaderDescriptor; bool m_HeaderDescriptorSet; }; } //end namespace mitk #endif // DIFFUSIONIMAGECREATIONFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.cxx b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.cxx similarity index 84% rename from Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.cxx rename to Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.cxx index 0779f47f46..7db96fb2e4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.cxx +++ b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.cxx @@ -1,244 +1,267 @@ /*=================================================================== 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 DIFFUSIONIMAGETRANSFORMEDCREATIONFILTER_CXX #define DIFFUSIONIMAGETRANSFORMEDCREATIONFILTER_CXX #include "mitkDiffusionImageTransformedCreationFilter.h" #include "mitkDiffusionImageCorrectionFilter.h" +#include "mitkImageCast.h" +#include "mitkImageWriteAccessor.h" +#include "mitkImageTimeSelector.h" +#include "mitkProperties.h" +#include "mitkIOUtil.h" + + #include #include #include +#include #include template static void ResampleImage( typename ItkImageType::Pointer itk_reference, mitk::Image::Pointer mitk_input, typename TTransformType::Pointer transform, unsigned int interpolator, mitk::Image::Pointer output_target, unsigned int position ) { typedef itk::LinearInterpolateImageFunction< ItkImageType, double > InterpolatorType; typename InterpolatorType::Pointer linear_interpolator = InterpolatorType::New(); typedef itk::NearestNeighborInterpolateImageFunction< ItkImageType, double > NearestNeighborInterpolatorType; typename NearestNeighborInterpolatorType::Pointer nn_interpolator = NearestNeighborInterpolatorType::New(); typedef itk::WindowedSincInterpolateImageFunction< ItkImageType, 7> WindowedSincInterpolatorType; typename WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New(); + + typedef itk::BSplineInterpolateImageFunction< ItkImageType, double, double> BSplineInterpolatorType; + typename BSplineInterpolatorType::Pointer bicubic_interpolator = BSplineInterpolatorType::New(); + bicubic_interpolator->SetSplineOrder(3); + typedef typename itk::ResampleImageFilter< ItkImageType, ItkImageType, double> ResampleImageFilterType; typename ResampleImageFilterType::Pointer resampler = ResampleImageFilterType::New(); - resampler->SetInterpolator( linear_interpolator ); - if( interpolator == mitk::DiffusionImageTransformedCreationFilter::WindowedSinc ) - resampler->SetInterpolator( sinc_interpolator ); - if ( interpolator == mitk::DiffusionImageTransformedCreationFilter::NearestNeighbor ) - resampler->SetInterpolator ( nn_interpolator ); + // Create interpolator pool + typedef itk::InterpolateImageFunction< ItkImageType, double> InterpolateFunctionBaseType; + std::vector< InterpolateFunctionBaseType* > InterpolateFunctionsPool; + + InterpolateFunctionsPool.push_back( nn_interpolator.GetPointer() ); + InterpolateFunctionsPool.push_back( linear_interpolator.GetPointer() ); + InterpolateFunctionsPool.push_back( bicubic_interpolator.GetPointer() ); + InterpolateFunctionsPool.push_back( sinc_interpolator.GetPointer() ); + + // select interpolator by the selection flag + resampler->SetInterpolator( InterpolateFunctionsPool.at( interpolator ) ); typename ItkImageType::Pointer itk_input; mitk::CastToItkImage< ItkImageType >( mitk_input, itk_input ); resampler->SetInput( itk_input ); resampler->SetTransform( transform ); resampler->SetReferenceImage( itk_reference ); resampler->UseReferenceImageOn(); resampler->Update(); mitk::Image::Pointer current_resampled = mitk::Image::New(); current_resampled->InitializeByItk( resampler->GetOutput() ); current_resampled->SetImportChannel( resampler->GetOutput()->GetBufferPointer(), mitk::Image::CopyMemory ); mitk::ImageWriteAccessor imac( current_resampled ); output_target->SetImportVolume( imac.GetData(), position, 0, mitk::Image::CopyMemory ); } /* template< typename TTransformType > static mitk::DiffusionPropertyHelper::GradientDirectionType TransformGradientDirection( mitk::DiffusionPropertyHelper::GradientDirectionType vec_in, typename TTransformType::Pointer transform ) { mitk::DiffusionPropertyHelper::GradientDirectionType vec_out; vec_out.fill(1.0); typedef typename TTransformType::MatrixType TMatrixType; return vec_out; } */ template< typename TTransformType> mitk::DiffusionImageTransformedCreationFilter ::DiffusionImageTransformedCreationFilter() { - + this->m_OutputPrefix = "/tmp/"; } template< typename TTransformType> mitk::DiffusionImageTransformedCreationFilter ::~DiffusionImageTransformedCreationFilter() { } template< typename TTransformType> mitk::DiffusionImageHeaderDescriptor mitk::DiffusionImageTransformedCreationFilter ::GetTransformedHeaderInformation() { typedef mitk::DiffusionPropertyHelper DPH; DiffusionImageHeaderDescriptor dwhdesc; dwhdesc.m_BValue = DPH::GetReferenceBValue( this->m_DiffusionReferenceImage ); // TODO : here comes transformation of the gradients dwhdesc.m_GradientDirections = DPH::GetGradientContainer( this->m_DiffusionReferenceImage ); dwhdesc.m_MeasurementFrame = DPH::GetMeasurementFrame( this->m_DiffusionReferenceImage ); dwhdesc.m_BValueMapType = DPH::GetBValueMap( this->m_DiffusionReferenceImage ); return dwhdesc; } +template< typename TTransformType> +void mitk::DiffusionImageTransformedCreationFilter< TTransformType > +::GenerateOutputInformation() +{ + MITK_INFO << "Debug info"; +} + + template< typename TTransformType> void mitk::DiffusionImageTransformedCreationFilter ::GenerateData() { mitk::Image::Pointer input = this->GetInput(); // validity checks if( m_InternalTransforms.size() != input->GetTimeSteps() ) { mitkThrow() << "Number of transforms" << m_InternalTransforms.size() << "differ from number of timesteps" << input->GetTimeSteps(); } typedef itk::Image< DiffusionPropertyHelper::DiffusionPixelType, 3> ImageType; // create intermediate output mitk::Image::Pointer resampled_output = mitk::Image::New(); resampled_output->Initialize( input ); ImageType::Pointer current_itk_reference = ImageType::New(); CastToItkImage( this->m_ResamplingReferenceImage, current_itk_reference ); - mitk::IOUtil::Save( input, "/tmp/resampled_3dnt_pre.nrrd"); - - auto time_dim = input->GetDimension(3); - for( auto time_idx = 0; time_idx < time_dim; time_idx++ ) + unsigned int time_dim = input->GetDimension(3); + for( unsigned int time_idx = 0; time_idx < time_dim; time_idx++ ) { mitk::ImageTimeSelector::Pointer t_sel = mitk::ImageTimeSelector::New(); t_sel->SetInput( input ); t_sel->SetTimeNr( time_idx ); t_sel->Update(); ResampleImage< TTransformType, ImageType>( current_itk_reference, t_sel->GetOutput(), this->m_InternalTransforms.at(time_idx), this->m_InterpolationLevel, resampled_output, time_idx); } // call here creation filter mitk::DiffusionImageCreationFilter::Pointer creator = mitk::DiffusionImageCreationFilter::New(); - mitk::IOUtil::Save( resampled_output, "/tmp/resampled_3dnt_post.nrrd"); - creator->SetInput( resampled_output ); creator->SetHeaderDescriptor( this->GetTransformedHeaderInformation() ); creator->Update(); - mitk::Image::Pointer output = this->GetOutput(); - output->Initialize( creator->GetOutput() ); + mitk::Image::Pointer output = creator->GetOutput(); typedef mitk::DiffusionPropertyHelper DPH; float BValue = mitk::DiffusionPropertyHelper::GetReferenceBValue( creator->GetOutput() ); mitk::BValueMapProperty::BValueMap BValueMap = mitk::BValueMapProperty::CreateBValueMap( mitk::DiffusionPropertyHelper::GetGradientContainer( creator->GetOutput() ), BValue ); output->SetProperty( DPH::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( mitk::DiffusionPropertyHelper::GetGradientContainer( creator->GetOutput() )) ); output->SetProperty( DPH::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( mitk::DiffusionPropertyHelper::GetMeasurementFrame( creator->GetOutput() ) ) ); output->SetProperty( DPH::BVALUEMAPPROPERTYNAME.c_str(), mitk::BValueMapProperty::New( BValueMap ) ); output->SetProperty( DPH::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( BValue ) ); // correct gradients mitk::DiffusionImageCorrectionFilter::Pointer corrector = mitk::DiffusionImageCorrectionFilter::New(); corrector->SetImage( output) ; corrector->CorrectDirections( this->m_RotationMatrices ); + DPH pHelper( output ); + pHelper.InitializeImage(); - output->Modified(); - - - mitk::IOUtil::Save( output, "/tmp/resampled_dw.dwi"); + m_OutputCache = output; + this->SetPrimaryOutput( m_OutputCache ); + m_OutputCache->Modified(); } template< typename TTransformType> void mitk::DiffusionImageTransformedCreationFilter ::SetTransformParameters(const TransformParametersContainerType ¶ms) { if( params.empty() ) { MITK_ERROR << "Empty parameter list given."; return; } TransformContainerType transforms; auto iter = std::begin( params ); while( iter != std::end( params ) ) { typename TTransformType::Pointer transform = TTransformType::New(); transform->SetParameters( (*iter) ); transforms.push_back( transform ); ++iter; } this->SetTransforms( transforms ); } template< typename TTransformType> void mitk::DiffusionImageTransformedCreationFilter ::SetTransforms(const TransformContainerType &transforms) { if( transforms.empty() ) { MITK_ERROR << "Empty transform list given."; return; } this->m_InternalTransforms.reserve( transforms.size() ); std::copy( transforms.begin(), transforms.end(), std::back_inserter(this->m_InternalTransforms ) ); MatrixType E; E.set_identity(); for( auto iter = std::begin(this->m_InternalTransforms); iter != std::end( this->m_InternalTransforms); ++iter) { MatrixType A = E * (*iter)->GetMatrix().GetVnlMatrix(); this->m_RotationMatrices.push_back( A ); } } #endif // DIFFUSIONIMAGETRANSFORMEDCREATIONFILTER_CXX diff --git a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.h index 93bbb0fd03..5039c443ce 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.h @@ -1,102 +1,116 @@ /*=================================================================== 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 DIFFUSIONIMAGETRANSFORMEDCREATIONFILTER_H #define DIFFUSIONIMAGETRANSFORMEDCREATIONFILTER_H #include "mitkDiffusionPropertyHelper.h" #include "mitkDiffusionImageCreationFilter.h" #include "mitkDiffusionImageCorrectionFilter.h" #include namespace mitk { /** * @brief DiffusionImageTransformedCreationFilter extends its superclass ( DiffusionImageCreationFilter ) by providing the possibility * to specify a list of transforms to be applied to the separate volumes before */ template class MITKDIFFUSIONCORE_EXPORT DiffusionImageTransformedCreationFilter : public ImageToImageFilter { public: /** class macros */ mitkClassMacro( DiffusionImageTransformedCreationFilter, ImageToImageFilter ) itkNewMacro(Self) enum InterpolatorLevel { NearestNeighbor = 0, - Linear, + Trilinear, + BSpline_Bicubic, WindowedSinc }; typedef std::vector< typename TTransformType::Pointer > TransformContainerType; typedef std::vector< typename TTransformType::ParametersType > TransformParametersContainerType; /** Set a transform container to be used, one transform per volume */ void SetTransforms( const TransformContainerType& ); /** Set a container of transform parameters to be used to create transforms for transforming the input */ void SetTransformParameters( const TransformParametersContainerType& ); /** Set interpolation type for resampling */ void SetInterpolationLevel( InterpolatorLevel lvl) { m_InterpolationLevel = lvl; } void GenerateData() override; void SetResamplingReferenceImage( mitk::Image::Pointer image ) { this->m_ResamplingReferenceImage = image; } void SetReferenceDiffusionImage( mitk::Image::Pointer dwimage ) { this->m_DiffusionReferenceImage = dwimage; } + void SetOutputPrefix( std::string prefix ) +{ + this->m_OutputPrefix = prefix; +} + + virtual void GenerateOutputInformation() override; + protected: DiffusionImageTransformedCreationFilter(); ~DiffusionImageTransformedCreationFilter(); DiffusionImageHeaderDescriptor GetTransformedHeaderInformation(); TransformContainerType m_InternalTransforms; mitk::Image::Pointer m_ResamplingReferenceImage; mitk::Image::Pointer m_DiffusionReferenceImage; InterpolatorLevel m_InterpolationLevel; typedef DiffusionImageCorrectionFilter::TransformMatrixType MatrixType; std::vector< MatrixType > m_RotationMatrices; -}; -#include "mitkDiffusionImageTransformedCreationFilter.cxx" + std::string m_OutputPrefix; + + mitk::Image::Pointer m_OutputCache; +}; } // end namespace mitk +#ifndef ITK_MANUAL_INSTANTIATION + #include "mitkDiffusionImageTransformedCreationFilter.cxx" +#endif + #endif // DIFFUSIONIMAGETRANSFORMEDCREATIONFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp index d93af7cf2d..e876e0d632 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp @@ -1,118 +1,117 @@ /*=================================================================== 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 MITKDIFFUSIONIMAGECORRECTIONFILTER_CPP #define MITKDIFFUSIONIMAGECORRECTIONFILTER_CPP #include "mitkDiffusionImageCorrectionFilter.h" #include #include #include mitk::DiffusionImageCorrectionFilter::DiffusionImageCorrectionFilter() { } mitk::DiffusionImageCorrectionFilter::TransformMatrixType mitk::DiffusionImageCorrectionFilter ::GetRotationComponent(const TransformMatrixType &A) { TransformMatrixType B; B = A * A.transpose(); // get the eigenvalues and eigenvectors typedef double MType; vnl_vector< MType > eigvals; vnl_matrix< MType > eigvecs; vnl_symmetric_eigensystem_compute< MType > ( B, eigvecs, eigvals ); vnl_matrix_fixed< MType, 3, 3 > eigvecs_fixed; eigvecs_fixed.set_columns(0, eigvecs ); TransformMatrixType C; C.set_identity(); vnl_vector_fixed< MType, 3 > eigvals_sqrt; for(unsigned int i=0; i<3; i++) { C(i,i) = std::sqrt( eigvals[i] ); } TransformMatrixType S = vnl_inverse( eigvecs_fixed * C * vnl_inverse( eigvecs_fixed )) * A; return S; } void mitk::DiffusionImageCorrectionFilter ::CorrectDirections( const TransformsVectorType& transformations) { if( m_SourceImage.IsNull() ) { mitkThrow() << " No diffusion image given! "; } GradientDirectionContainerPointerType directions = static_cast( m_SourceImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); GradientDirectionContainerPointerType corrected_directions = GradientDirectionContainerType::New(); unsigned int transformed = 0; for(size_t i=0; i< directions->Size(); i++ ) { // skip b-zero images if( directions->ElementAt(i).one_norm() <= 0.0 ) { corrected_directions->push_back( directions->ElementAt(i) ); continue; } GradientDirectionType corrected = GetRotationComponent( transformations.at(transformed)) * directions->ElementAt(i); // store the corrected direction corrected_directions->push_back( corrected ); transformed++; } // replace the old directions with the corrected ones m_SourceImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( corrected_directions ) ); - } void mitk::DiffusionImageCorrectionFilter ::CorrectDirections( const TransformMatrixType& transformation) { if( m_SourceImage.IsNull() ) { mitkThrow() << " No diffusion image given! "; } TransformsVectorType transfVec; for (unsigned int i=0; i< static_cast( m_SourceImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->Size();i++) { transfVec.push_back(transformation); } this->CorrectDirections(transfVec); } #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.cpp index 828cd1a7c4..a5154d5ddd 100644 --- a/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/src/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.cpp @@ -1,212 +1,231 @@ /*=================================================================== 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 "mitkDiffusionImageCreationFilter.h" #include "mitkProperties.h" #include "mitkImageTimeSelector.h" #include "mitkImageCast.h" #include "mitkImageToItk.h" #include "mitkImageAccessByItk.h" #include "mitkITKImageImport.h" #include "mitkIOUtil.h" +#include + #include /** * @brief RemapIntoVectorImage Take a 3d+t image and reinterpret it as vector image * @return vectoriamge */ -mitk::DiffusionImageCreationFilter::VectorImageType::Pointer +mitk::Image::Pointer mitk::DiffusionImageCreationFilter::RemapIntoVectorImage( mitk::Image::Pointer input) { typedef itk::Image ImageVolumeType; typedef itk::ComposeImageFilter< ImageVolumeType > ComposeFilterType; ComposeFilterType::Pointer vec_composer = ComposeFilterType::New(); mitk::ImageTimeSelector::Pointer t_selector = mitk::ImageTimeSelector::New(); t_selector->SetInput( input ); for( unsigned int i=0; i< input->GetTimeSteps(); i++) { t_selector->SetTimeNr(i); t_selector->Update(); ImageVolumeType::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(); + MITK_ERROR << "Caught exception while updating compose filter: " << e.what(); } mitk::DiffusionImageCreationFilter::VectorImageType::Pointer vector_image = vec_composer->GetOutput(); - vector_image->GetPixelContainer()->SetContainerManageMemory(false); + m_OutputCache = mitk::Image::New(); + + //mitk::GrabItkImageMemory( vector_image ); + m_OutputCache->InitializeByItk( vector_image.GetPointer() ); + m_OutputCache->SetImportVolume( vector_image->GetBufferPointer(), 0, 0, Image::CopyMemory ); + vector_image->GetPixelContainer()->ContainerManageMemoryOff(); - return vector_image; + return m_OutputCache; } mitk::DiffusionImageCreationFilter::DiffusionImageCreationFilter() : m_ReferenceImage( nullptr ) { m_HeaderDescriptorSet = false; this->SetNumberOfRequiredInputs(1); this->SetNumberOfRequiredOutputs(1); } mitk::DiffusionImageCreationFilter::~DiffusionImageCreationFilter() { } +void mitk::DiffusionImageCreationFilter +::GenerateOutputInformation() +{ + mitk::Image::Pointer input = this->GetInput(); + mitk::Image::Pointer output = this->GetOutput(); + + // the filter merges all 3d+t to the components -> we need to adapt the pixel type accordingly + output->Initialize( MakePixelType< mitk::DiffusionPropertyHelper::DiffusionPixelType, 3>( input->GetTimeSteps() ), + *input->GetGeometry(0) ); +} + void mitk::DiffusionImageCreationFilter::SetReferenceImage( mitk::Image::Pointer reference_image ) { if( reference_image.IsNull() ) { mitkThrow() << "Null-pointer image provided as reference. "; } if( ! DPH::IsDiffusionWeightedImage(reference_image) ) { mitkThrow() << "The image provided as reference is not a diffusion-weighted image. Cannot proceed. "; } this->m_ReferenceImage = reference_image; } void mitk::DiffusionImageCreationFilter::GenerateData() { const mitk::Image::Pointer input_image = this->GetInput(0); if( input_image.IsNull() ) { mitkThrow() << "No input specified. Cannot proceed "; } if( !( m_HeaderDescriptorSet ^ m_ReferenceImage.IsNotNull() ) ) { mitkThrow() << "Either a header descriptor or a reference diffusion-weighted image have to be provided. Terminating."; } - // data packed into vector image - OutputType::Pointer outputForCache = this->GetOutput(); + mitk::Image::Pointer outputForCache = RemapIntoVectorImage( input_image ); + - mitk::Image::Pointer mitkvectorimage = mitk::GrabItkImageMemory( RemapIntoVectorImage( input_image )); - outputForCache->Initialize( mitkvectorimage ); // header information GradientDirectionContainerType::Pointer DiffusionVectors = this->InternalGetGradientDirections( ); MeasurementFrameType MeasurementFrame = this->InternalGetMeasurementFrame(); float BValue = this->InternalGetBValue(); // create BValueMap mitk::BValueMapProperty::BValueMap BValueMap = mitk::BValueMapProperty::CreateBValueMap( DiffusionVectors, BValue ); - outputForCache->SetProperty( DPH::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( DiffusionVectors ) ); - outputForCache->SetProperty( DPH::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( MeasurementFrame ) ); - outputForCache->SetProperty( DPH::BVALUEMAPPROPERTYNAME.c_str(), mitk::BValueMapProperty::New( BValueMap ) ); - outputForCache->SetProperty( DPH::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( BValue ) ); + m_OutputCache->SetProperty( DPH::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( DiffusionVectors ) ); + m_OutputCache->SetProperty( DPH::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( MeasurementFrame ) ); + m_OutputCache->SetProperty( DPH::BVALUEMAPPROPERTYNAME.c_str(), mitk::BValueMapProperty::New( BValueMap ) ); + m_OutputCache->SetProperty( DPH::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( BValue ) ); + + DPH pHelper( m_OutputCache ); + pHelper.InitializeImage(); - outputForCache->Modified(); + this->SetPrimaryOutput( m_OutputCache ); } void mitk::DiffusionImageCreationFilter::SetHeaderDescriptor(DiffusionImageHeaderDescriptor header_descriptor) { this->m_HeaderDescriptor = header_descriptor; this->m_HeaderDescriptorSet = true; } mitk::DiffusionImageCreationFilter::MeasurementFrameType mitk::DiffusionImageCreationFilter::InternalGetMeasurementFrame() { MeasurementFrameType MeasurementFrame; if( m_ReferenceImage.IsNotNull() ) { MeasurementFrame = DPH::GetMeasurementFrame( m_ReferenceImage ); } else if ( m_HeaderDescriptorSet ) { MeasurementFrame = m_HeaderDescriptor.m_MeasurementFrame; } else { MeasurementFrame(0,0) = 1; MeasurementFrame(0,1) = 0; MeasurementFrame(0,2) = 0; MeasurementFrame(1,0) = 0; MeasurementFrame(1,1) = 1; MeasurementFrame(1,2) = 0; MeasurementFrame(2,0) = 0; MeasurementFrame(2,1) = 0; MeasurementFrame(2,2) = 1; MITK_WARN << "Created default measurement frame as non provided ( no reference image or header information provided)"; } return MeasurementFrame; } mitk::DiffusionImageCreationFilter::GradientDirectionContainerType::Pointer mitk::DiffusionImageCreationFilter::InternalGetGradientDirections() { GradientDirectionContainerType::Pointer DiffusionVectors = GradientDirectionContainerType::New(); if( this->m_ReferenceImage ) { DiffusionVectors = DPH::GetGradientContainer( this->m_ReferenceImage ); } else if ( m_HeaderDescriptorSet ) { DiffusionVectors = m_HeaderDescriptor.m_GradientDirections; } return DiffusionVectors; } float mitk::DiffusionImageCreationFilter::InternalGetBValue() { float bvalue = -1; if( m_ReferenceImage.IsNotNull() ) { bvalue = DPH::GetReferenceBValue( m_ReferenceImage ); } else if ( m_HeaderDescriptorSet ) { bvalue = m_HeaderDescriptor.m_BValue; } else { MITK_ERROR << "No reference image and no header descriptor provided."; } return bvalue; }