diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h index 6e7bccd182..f1e55fa6ac 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h @@ -1,220 +1,220 @@ /*=================================================================== 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 __itkDiffusionMultiShellQballReconstructionImageFilter_h_ #define __itkDiffusionMultiShellQballReconstructionImageFilter_h_ #include namespace itk{ /** \class DiffusionMultiShellQballReconstructionImageFilter I. Aganj, C. Lenglet, G. Sapiro, E. Yacoub, K. Ugurbil, and N. Harel, “Reconstruction of the orientation distribution function in single and multiple shell q-ball imaging within constant solid angle,” Magnetic Resonance in Medicine, vol. 64, no. 2, pp. 554–566, 2010. */ template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> class DiffusionMultiShellQballReconstructionImageFilter : public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > { public: typedef DiffusionMultiShellQballReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef TReferenceImagePixelType ReferencePixelType; /** GradientImageType * (e.g. type short)*/ typedef TGradientImagePixelType GradientPixelType; /** GradientImageType * 3D VectorImage containing GradientPixelTypes */ typedef VectorImage< GradientPixelType, 3 > GradientImagesType; /** ODF PixelType */ typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType; /** ODF ImageType */ typedef Image< OdfPixelType, 3 > OdfImageType; /** BzeroImageType */ typedef Image< TOdfPixelType, 3 > BZeroImageType; /** Container to hold gradient directions of the 'n' DW measurements */ typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType; typedef Image< Vector< TOdfPixelType, (NOrderL*NOrderL + NOrderL + 2)/2 + NOrderL >, 3 > CoefficientImageType; typedef std::map > BValueMap; typedef std::map >::iterator BValueMapIteraotr; typedef std::vector IndiciesVector; // --------------------------------------------------------------------------------------------// /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(DiffusionMultiShellQballReconstructionImageFilter, ImageToImageFilter) /** Get reference image */ virtual typename Superclass::InputImageType * GetInputImage() { return ( static_cast< typename Superclass::InputImageType *>(this->ProcessObject::GetInput(0)) ); } /** Replaces the Input method. - * Var vols = mitk::DiffusionImage + * Var vols = mitk::Image * ----------------------------------------------------- - * GradientDirectionContainerType-Input gradientDirectionContainer (e.g. vols->GetDirections) - * GradientImagesType-Input gradientImage (e.g. vols->GetVectorImage) - * float-Input bvalue (e.g. vols->GetB_Value) */ + * GradientDirectionContainerType-Input gradientDirectionContainer (e.g. GradientDirectionsContainerProperty) + * GradientImagesType-Input gradientImage (e.g. itkVectorImage) + * float-Input bvalue (e.g. ReferenceBValueProperty) */ void SetGradientImage( const GradientDirectionContainerType * gradientDirectionContainer, const GradientImagesType *gradientImage , float bvalue);//, std::vector listOfUserSelctedBValues ); /** Set a BValue Map (key = bvalue, value = indicies splittet for each shell) * If the input image containes more than three q-shells * (e.g. b-Values of 0, 1000, 2000, 3000, 4000, ...). * For the Analytical-Reconstruction it is needed to set a * BValue Map containing three shells in an arithmetic series * (e.g. 0, 1000, 2000, 3000). */ inline void SetBValueMap(BValueMap map){this->m_BValueMap = map;} /** Threshold on the reference image data. The output ODF will be a null * pdf for pixels in the reference image that have a value less than this * threshold. */ itkSetMacro( Threshold, ReferencePixelType ) itkGetMacro( Threshold, ReferencePixelType ) itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer ) /** Return non-diffusion weighted images */ itkGetMacro( BZeroImage, typename BZeroImageType::Pointer) /** Factor for Laplacian-Baltrami smoothing of the SH-coefficients*/ itkSetMacro( Lambda, double ) itkGetMacro( Lambda, double ) protected: DiffusionMultiShellQballReconstructionImageFilter(); ~DiffusionMultiShellQballReconstructionImageFilter() { } void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads ); private: enum ReconstructionType { Mode_Analytical3Shells, Mode_NumericalNShells, Mode_Standard1Shell }; ReconstructionType m_ReconstructionType; // Interpolation bool m_Interpolation_Flag; vnl_matrix< double > * m_Interpolation_SHT1_inv; vnl_matrix< double > * m_Interpolation_SHT2_inv; vnl_matrix< double > * m_Interpolation_SHT3_inv; vnl_matrix< double > * m_TARGET_SH_shell1; vnl_matrix< double > * m_TARGET_SH_shell2; vnl_matrix< double > * m_TARGET_SH_shell3; unsigned int m_MaxDirections; vnl_matrix< double > * m_CoeffReconstructionMatrix; vnl_matrix< double > * m_ODFSphericalHarmonicBasisMatrix; /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; /** Number of baseline images */ unsigned int m_NumberOfBaselineImages; /** Threshold on the reference image data */ ReferencePixelType m_Threshold; typename BZeroImageType::Pointer m_BZeroImage; typename CoefficientImageType::Pointer m_CoefficientImage; float m_BValue; BValueMap m_BValueMap; double m_Lambda; bool m_IsHemisphericalArrangementOfGradientDirections; bool m_IsArithmeticProgession; void ComputeReconstructionMatrix(IndiciesVector const & refVector); void ComputeODFSHBasis(); bool CheckDuplicateDiffusionGradients(); bool CheckForDifferingShellDirections(); IndiciesVector GetAllDirections(); void ComputeSphericalHarmonicsBasis(vnl_matrix* QBallReference, vnl_matrix* SHBasisOutput, int Lorder , vnl_matrix* LaplaciaBaltramiOutput =0 , vnl_vector* SHOrderAssociation =0 , vnl_matrix * SHEigenvalues =0); void Normalize(OdfPixelType & odf ); void S_S0Normalization( vnl_vector & vec, double b0 = 0 ); void DoubleLogarithm(vnl_vector & vec); double CalculateThreashold(const double value, const double delta); void Projection1(vnl_vector & vec, double delta = 0.01); void Projection2( vnl_vector & E1, vnl_vector & E2, vnl_vector & E3, double delta = 0.01); void Projection3( vnl_vector & A, vnl_vector & alpha, vnl_vector & beta, double delta = 0.01); void StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread); void AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread); void NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread); void GenerateAveragedBZeroImage(const OutputImageRegionType& outputRegionForThread); void ComputeSphericalFromCartesian(vnl_matrix * Q, const IndiciesVector & refShell); //------------------------- VNL-function ------------------------------------ template vnl_vector< WntValue> element_cast (vnl_vector< CurrentValue> const& v1) { vnl_vector result(v1.size()); for(unsigned int i = 0 ; i < v1.size(); i++) result[i] = static_cast< WntValue>(v1[i]); return result; } template double dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ) { double result = (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (v1.two_norm() * v2.two_norm()); return result ; } }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h index 32bccf0800..f465a7c29c 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h @@ -1,117 +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 __mitkTeemDiffusionTensor3DReconstructionImageFilter_h__ #define __mitkTeemDiffusionTensor3DReconstructionImageFilter_h__ #include "mitkImage.h" #include "mitkTensorImage.h" -#include "mitkDiffusionImage.h" #include "itkDiffusionTensor3D.h" namespace mitk { enum TeemTensorEstimationMethods{ TeemTensorEstimationMethodsLLS, TeemTensorEstimationMethodsNLS, TeemTensorEstimationMethodsWLS, TeemTensorEstimationMethodsMLE, }; template< class DiffusionImagePixelType = short, class TTensorPixelType=float > class TeemDiffusionTensor3DReconstructionImageFilter : public itk::Object { public: typedef TTensorPixelType TensorPixelType; typedef itk::Vector VectorType; typedef itk::Image VectorImageType; typedef itk::DiffusionTensor3D TensorType; typedef itk::Image ItkTensorImageType; typedef itk::Vector ItkTensorVectorType; typedef itk::Image ItkTensorVectorImageType; typedef DiffusionImagePixelType DiffusionPixelType; typedef itk::VectorImage< DiffusionPixelType, 3 > DiffusionImageType; mitkClassMacro( TeemDiffusionTensor3DReconstructionImageFilter, itk::Object ); itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkGetMacro(Input, - typename DiffusionImage::Pointer); + mitk::Image::Pointer); itkSetMacro(Input, - typename DiffusionImage::Pointer); + mitk::Image::Pointer); itkGetMacro(EstimateErrorImage, bool); itkSetMacro(EstimateErrorImage, bool); itkGetMacro(Sigma, float); itkSetMacro(Sigma, float); itkGetMacro(EstimationMethod, TeemTensorEstimationMethods); itkSetMacro(EstimationMethod, TeemTensorEstimationMethods); itkGetMacro(NumIterations, int); itkSetMacro(NumIterations, int); itkGetMacro(ConfidenceThreshold, double); itkSetMacro(ConfidenceThreshold, double); itkGetMacro(ConfidenceFuzzyness, float); itkSetMacro(ConfidenceFuzzyness, float); itkGetMacro(MinPlausibleValue, double); itkSetMacro(MinPlausibleValue, double); itkGetMacro(Output, mitk::TensorImage::Pointer); itkGetMacro(OutputItk, mitk::TensorImage::Pointer); // do the work virtual void Update(); protected: TeemDiffusionTensor3DReconstructionImageFilter(); virtual ~TeemDiffusionTensor3DReconstructionImageFilter(); - typename DiffusionImage::Pointer m_Input; + mitk::Image::Pointer m_Input; bool m_EstimateErrorImage; float m_Sigma; TeemTensorEstimationMethods m_EstimationMethod; int m_NumIterations; double m_ConfidenceThreshold; float m_ConfidenceFuzzyness; double m_MinPlausibleValue; mitk::TensorImage::Pointer m_Output; mitk::TensorImage::Pointer m_OutputItk; mitk::Image::Pointer m_ErrorImage; }; } #include "mitkTeemDiffusionTensor3DReconstructionImageFilter.cpp" #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp index d11c505925..3f6fbc61a9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp @@ -1,327 +1,316 @@ /*=================================================================== 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 +#ifndef MITKDWIHEADMOTIONCORRECTIONFILTER_CPP +#define MITKDWIHEADMOTIONCORRECTIONFILTER_CPP #include "mitkDWIHeadMotionCorrectionFilter.h" #include "itkSplitDWImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "mitkImageTimeSelector.h" #include "mitkPyramidImageRegistrationMethod.h" //#include "mitkRegistrationMethodITK4.h" -#include "mitkImageToDiffusionImageSource.h" +#include #include "mitkDiffusionImageCorrectionFilter.h" #include +#include #include #include "mitkIOUtil.h" #include -template< typename DiffusionPixelType> -mitk::DWIHeadMotionCorrectionFilter - ::DWIHeadMotionCorrectionFilter() +mitk::DWIHeadMotionCorrectionFilter::DWIHeadMotionCorrectionFilter() : m_CurrentStep(0), m_Steps(100), m_IsInValidState(true), m_AbortRegistration(false), m_AverageUnweighted(true) { } -template< typename DiffusionPixelType> -void mitk::DWIHeadMotionCorrectionFilter -::GenerateData() + +void mitk::DWIHeadMotionCorrectionFilter::GenerateData() { typedef itk::SplitDWImageFilter< DiffusionPixelType, DiffusionPixelType> SplitFilterType; - DiffusionImageType* input = const_cast(this->GetInput(0)); + InputImageType* input = const_cast(this->GetInput(0)); + + ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); + mitk::CastToItkImage(input, itkVectorImagePointer); typedef mitk::PyramidImageRegistrationMethod RegistrationMethod; // typedef mitk::RegistrationMethodITKV4 RegistrationMethod - - unsigned int numberOfSteps = input->GetVectorImage()->GetNumberOfComponentsPerPixel () ; + unsigned int numberOfSteps = itkVectorImagePointer->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() ); + B0ExtractorType::Pointer b0_extractor = B0ExtractorType::New(); + b0_extractor->SetInput( itkVectorImagePointer ); + b0_extractor->SetDirections( static_cast( input->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ); b0_extractor->Update(); mitk::Image::Pointer b0Image = mitk::Image::New(); b0Image->InitializeByItk( b0_extractor->GetOutput() ); b0Image->SetImportChannel( b0_extractor->GetOutput()->GetBufferPointer(), mitk::Image::CopyMemory ); // (2.1) Use the extractor to access the extracted b0 volumes mitk::ImageTimeSelector::Pointer t_selector = 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(); // register b-zeros only if the flag to average is set, use the first one otherwise if( m_AverageUnweighted && numberOfb0Images > 1) { RegistrationMethod::Pointer registrationMethod = RegistrationMethod::New(); registrationMethod->SetFixedImage( b0referenceImage ); registrationMethod->SetTransformToRigid(); // the unweighted images are of same modality registrationMethod->SetCrossModalityOff(); // use the advanced (windowed sinc) interpolation registrationMethod->SetUseAdvancedInterpolation(true); // Initialize the temporary output image mitk::Image::Pointer registeredB0Image = b0Image->Clone(); mitk::ImageTimeSelector::Pointer t_selector2 = mitk::ImageTimeSelector::New(); t_selector2->SetInput( b0Image ); for( unsigned int i=1; iSetTimeNr(i); t_selector2->Update(); 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(); } // import volume to the inter-results mitk::ImageWriteAccessor imac(registrationMethod->GetResampledMovingImage()); registeredB0Image->SetImportVolume( imac.GetData(), i, 0, mitk::Image::ReferenceMemory ); } // use the accumulateImageFilter as provided by the ItkAccumulateFilter method in the header file AccessFixedDimensionByItk_1(registeredB0Image, ItkAccumulateFilter, (4), b0referenceImage ); } // // (2) Split the diffusion image into a 3d+t regular image, extract only the weighted images // - typename SplitFilterType::Pointer split_filter = SplitFilterType::New(); - split_filter->SetInput (input->GetVectorImage() ); - split_filter->SetExtractAllAboveThreshold(8, input->GetBValueMap() ); + SplitFilterType::Pointer split_filter = SplitFilterType::New(); + split_filter->SetInput (itkVectorImagePointer ); + split_filter->SetExtractAllAboveThreshold(8, static_cast(input->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap() ); 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 ); // // (3) Use again the time-selector to access the components separately in order // to perform the registration of Image -> unweighted reference // 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); weightedRegistrationMethod->SetVerboseOn(); // // - (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 // in own scope to release the accessor asap after copy { mitk::ImageWriteAccessor imac(b0referenceImage); registeredWeighted->SetImportVolume( imac.GetData(), 0,0, mitk::Image::CopyMemory ); } // mitk::Image::Pointer registeredWeighted = splittedImage->Clone(); // this time start at 0, we have only gradient images in the 3d+t file // the reference image comes form an other image mitk::ImageTimeSelector::Pointer t_selector_w = mitk::ImageTimeSelector::New(); t_selector_w->SetInput( splittedImage ); // store the rotation parts of the transformations in a vector typedef RegistrationMethod::TransformMatrixType MatrixType; std::vector< MatrixType > estimated_transforms; for( unsigned int i=0; iSetTimeNr(i); t_selector_w->Update(); weightedRegistrationMethod->SetMovingImage( t_selector_w->GetOutput() ); try { MITK_INFO << " === (" << i+1 <<"/"<< maxImageIdx << ") :: Starting registration"; weightedRegistrationMethod->Update(); } catch( const itk::ExceptionObject& e) { m_IsInValidState = false; mitkThrow() << "Failed to register the b0 images, the PyramidRegistration threw an exception: \n" << e.what(); } // allow expansion mitk::ImageWriteAccessor imac(weightedRegistrationMethod->GetResampledMovingImage()); registeredWeighted->SetImportVolume( imac.GetData(), i+1, 0, mitk::Image::CopyMemory); estimated_transforms.push_back( weightedRegistrationMethod->GetLastRotationMatrix() ); } // // (4) Cast the resulting image back to an diffusion weighted image // - typename DiffusionImageType::GradientDirectionContainerType *gradients = input->GetDirections(); - typename DiffusionImageType::GradientDirectionContainerType::Pointer gradients_new = - DiffusionImageType::GradientDirectionContainerType::New(); - typename DiffusionImageType::GradientDirectionType bzero_vector; + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType *gradients = static_cast( input->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradients_new = + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New(); + mitk::DiffusionPropertyHelper::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(); + SplitFilterType::IndexListType index_list = split_filter->GetIndexList(); + 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(); - - caster->SetImage( registeredWeighted ); - caster->SetBValue( input->GetReferenceBValue() ); - caster->SetGradientDirections( gradients_new.GetPointer() ); - - try - { - caster->Update(); - } - catch( const itk::ExceptionObject& e) - { - m_IsInValidState = false; - MITK_ERROR << "Casting back to diffusion image failed: "; - mitkThrow() << "Subprocess failed with exception: " << e.what(); - } + registeredWeighted->SetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( gradients_new.GetPointer() )); + registeredWeighted->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(input->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); // // (5) Adapt the gradient directions according to the estimated transforms // - typedef mitk::DiffusionImageCorrectionFilter< DiffusionPixelType > CorrectionFilterType; - typename CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); + typedef mitk::DiffusionImageCorrectionFilter CorrectionFilterType; + CorrectionFilterType::Pointer corrector = CorrectionFilterType::New(); - OutputImagePointerType output = caster->GetOutput(); + mitk::Image::Pointer output = registeredWeighted; corrector->SetImage( output ); corrector->CorrectDirections( estimated_transforms ); // // (6) Pass the corrected image to the filters output port // m_CurrentStep += 1; - this->GetOutput()->SetVectorImage(output->GetVectorImage()); - this->GetOutput()->SetReferenceBValue(output->GetReferenceBValue()); - this->GetOutput()->SetMeasurementFrame(output->GetMeasurementFrame()); - this->GetOutput()->SetDirections(output->GetDirections()); - this->GetOutput()->InitializeFromVectorImage(); + + this->GetOutput()->Initialize( output ); + + this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( static_cast( output->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer() ) ); + this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( static_cast( output->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() )->GetMeasurementFrame() ) ); + this->GetOutput()->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(output->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); + + mitk::DiffusionPropertyHelper propertyHelper( this->GetOutput() ); + propertyHelper.InitializeImage(); this->GetOutput()->Modified(); } -#endif // MITKDIFFUSIONIMAGETODIFFUSIONIMAGEFILTER_CPP +#endif // MITKDWIHEADMOTIONCORRECTIONFILTER_CPP diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h index d8065ff27d..e159597b08 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h @@ -1,132 +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 "mitkImageToImageFilter.h" #include #include #include "mitkITKImageImport.h" +#include +#include 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 > + +class MitkDiffusionCore_EXPORT DWIHeadMotionCorrectionFilter + : public ImageToImageFilter { public: // class macros mitkClassMacro( DWIHeadMotionCorrectionFilter, - DiffusionImageToDiffusionImageFilter ) + ImageToImageFilter ) itkFactorylessNewMacro(Self) itkCloneMacro(Self) itkGetMacro( Steps, unsigned long ) itkGetMacro( CurrentStep, unsigned long ) itkGetMacro( IsInValidState, bool) itkSetMacro( AbortRegistration, bool ) itkSetMacro( AverageUnweighted, bool ) // public typedefs - typedef typename Superclass::InputImageType DiffusionImageType; - typedef typename Superclass::InputImagePointerType DiffusionImagePointerType; - - typedef typename Superclass::OutputImageType OutputImageType; - typedef typename Superclass::OutputImagePointerType OutputImagePointerType; + typedef short DiffusionPixelType; + typedef Superclass::InputImageType InputImageType; + typedef Superclass::OutputImageType OutputImageType; + typedef itk::VectorImage ITKDiffusionImageType; protected: DWIHeadMotionCorrectionFilter(); virtual ~DWIHeadMotionCorrectionFilter() {} virtual void GenerateData(); unsigned long m_CurrentStep; unsigned long m_Steps; bool m_IsInValidState; ///< Whether the filter is in a valid state, false if error occured bool m_AbortRegistration; ///< set flag to abort bool m_AverageUnweighted; }; /** * @brief Averages an 3d+t image along the time axis. * * The method uses the AccumulateImageFilter as provided by ITK and collapses the given 3d+t image * to an 3d image while computing the average over the time axis for each of the spatial voxels. */ template< typename TPixel, unsigned int VDimensions> static void ItkAccumulateFilter( const itk::Image< TPixel, VDimensions>* image, mitk::Image::Pointer& output) { // 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/mitkRegistrationWrapper.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkRegistrationWrapper.cpp index 0e1982705b..2f488f4334 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkRegistrationWrapper.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Registration/mitkRegistrationWrapper.cpp @@ -1,196 +1,216 @@ #include "mitkRegistrationWrapper.h" #include "mitkPyramidImageRegistrationMethod.h" -#include "mitkDiffusionImage.h" +#include "mitkImage.h" #include #include "itkB0ImageExtractionImageFilter.h" #include #include #include #include #include +#include #include void mitk::RegistrationWrapper::ApplyTransformationToImage(mitk::Image::Pointer img, const mitk::RegistrationWrapper::RidgidTransformType &transformation,double* offset, mitk::Image* resampleReference, bool binary) { - typedef mitk::DiffusionImage DiffusionImageType; + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientDirections = + static_cast( img->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); - if (dynamic_cast (img.GetPointer()) == NULL) + if (gradientDirections.IsNull() || ( gradientDirections->Size() == 0) ) { ItkImageType::Pointer itkImage = ItkImageType::New(); CastToItkImage(img, itkImage); typedef itk::Euler3DTransform< double > RigidTransformType; RigidTransformType::Pointer rtransform = RigidTransformType::New(); RigidTransformType::ParametersType parameters(RigidTransformType::ParametersDimension); for (int i = 0; i<6;++i) parameters[i] = transformation[i]; rtransform->SetParameters( parameters ); mitk::Point3D origin = itkImage->GetOrigin(); origin[0]-=offset[0]; origin[1]-=offset[1]; origin[2]-=offset[2]; mitk::Point3D newOrigin = rtransform->GetInverseTransform()->TransformPoint(origin); itk::Matrix dir = itkImage->GetDirection(); itk::Matrix transM ( vnl_inverse(rtransform->GetMatrix().GetVnlMatrix())); itk::Matrix newDirection = transM * dir; itkImage->SetOrigin(newOrigin); itkImage->SetDirection(newDirection); // Perform Resampling if reference image is provided if (resampleReference != NULL) { typedef itk::ResampleImageFilter ResampleFilterType; ItkImageType::Pointer itkReference = ItkImageType::New(); CastToItkImage(resampleReference,itkReference); typedef itk::Function::WelchWindowFunction<4> WelchWindowFunction; typedef itk::WindowedSincInterpolateImageFunction< ItkImageType, 4,WelchWindowFunction> WindowedSincInterpolatorType; WindowedSincInterpolatorType::Pointer sinc_interpolator = WindowedSincInterpolatorType::New(); typedef itk::LinearInterpolateImageFunction< ItkImageType> LinearInterpolatorType; LinearInterpolatorType::Pointer lin_interpolator = LinearInterpolatorType::New(); typedef itk::NearestNeighborInterpolateImageFunction< ItkImageType, double > NearestNeighborInterpolatorType; NearestNeighborInterpolatorType::Pointer nn_interpolator = NearestNeighborInterpolatorType::New(); typedef itk::BSplineInterpolateImageFunction< ItkImageType, double > BSplineInterpolatorType; BSplineInterpolatorType::Pointer bSpline_interpolator = BSplineInterpolatorType::New(); ResampleFilterType::Pointer resampler = ResampleFilterType::New(); resampler->SetInput(itkImage); resampler->SetReferenceImage( itkReference ); resampler->UseReferenceImageOn(); if (binary) resampler->SetInterpolator(nn_interpolator); else resampler->SetInterpolator(lin_interpolator); resampler->Update(); GrabItkImageMemory(resampler->GetOutput(), img); } else { // !! CastToItk behaves very differently depending on the original data type // if the target type is the same as the original, only a pointer to the data is set // and an additional GrabItkImageMemory will cause a segfault when the image is destroyed // GrabItkImageMemory - is not necessary in this case since we worked on the original data // See Bug 17538. if (img->GetPixelType().GetComponentTypeAsString() != "double") img = GrabItkImageMemory(itkImage); } } else { - DiffusionImageType::Pointer diffImages = dynamic_cast(img.GetPointer()); + mitk::Image::Pointer diffImages = dynamic_cast(img.GetPointer()); typedef itk::Euler3DTransform< double > RigidTransformType; RigidTransformType::Pointer rtransform = RigidTransformType::New(); RigidTransformType::ParametersType parameters(RigidTransformType::ParametersDimension); for (int i = 0; i<6;++i) { parameters[i] = transformation[i]; } rtransform->SetParameters( parameters ); - mitk::Point3D b0origin = diffImages->GetVectorImage()->GetOrigin(); + typedef itk::VectorImage ITKDiffusionImageType; + ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); + mitk::CastToItkImage(diffImages, itkVectorImagePointer); + + mitk::Point3D b0origin = itkVectorImagePointer->GetOrigin(); b0origin[0]-=offset[0]; b0origin[1]-=offset[1]; b0origin[2]-=offset[2]; mitk::Point3D newOrigin = rtransform->GetInverseTransform()->TransformPoint(b0origin); - itk::Matrix dir = diffImages->GetVectorImage()->GetDirection(); + itk::Matrix dir = itkVectorImagePointer->GetDirection(); itk::Matrix transM ( vnl_inverse(rtransform->GetMatrix().GetVnlMatrix())); itk::Matrix newDirection = transM * dir; - diffImages->GetVectorImage()->SetOrigin(newOrigin); - diffImages->GetVectorImage()->SetDirection(newDirection); + itkVectorImagePointer->SetOrigin(newOrigin); + itkVectorImagePointer->SetDirection(newDirection); diffImages->Modified(); - mitk::DiffusionImageCorrectionFilter::Pointer correctionFilter = mitk::DiffusionImageCorrectionFilter::New(); + mitk::DiffusionImageCorrectionFilter::Pointer correctionFilter = mitk::DiffusionImageCorrectionFilter::New(); // For Diff. Images: Need to rotate the gradients (works in-place) correctionFilter->SetImage(diffImages); correctionFilter->CorrectDirections(transM.GetVnlMatrix()); img = diffImages; } } void mitk::RegistrationWrapper::GetTransformation(mitk::Image::Pointer fixedImage, mitk::Image::Pointer movingImage, RidgidTransformType transformation,double* offset, bool useSameOrigin, mitk::Image* mask) { // Handle the case that fixed/moving image is a DWI image - mitk::DiffusionImage* fixedDwi = dynamic_cast*> (fixedImage.GetPointer()); - mitk::DiffusionImage* movingDwi = dynamic_cast*> (movingImage.GetPointer()); + mitk::Image* fixedDwi = dynamic_cast (fixedImage.GetPointer()); + mitk::Image* movingDwi = dynamic_cast (movingImage.GetPointer()); itk::B0ImageExtractionImageFilter::Pointer b0Extraction = itk::B0ImageExtractionImageFilter::New(); offset[0]=offset[1]=offset[2]=0; - if (fixedDwi != NULL) + + typedef itk::VectorImage ITKDiffusionImageType; + ITKDiffusionImageType::Pointer itkFixedDwiPointer = ITKDiffusionImageType::New(); + mitk::CastToItkImage(fixedDwi, itkFixedDwiPointer); + + ITKDiffusionImageType::Pointer itkMovingDwiPointer = ITKDiffusionImageType::New(); + mitk::CastToItkImage(movingDwi, itkMovingDwiPointer); + + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer fixedGradientDirections = + static_cast( fixedDwi->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); + + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer movingGradientDirections = + static_cast( movingDwi->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); + + if (fixedGradientDirections.IsNull() || ( fixedGradientDirections->Size() == 0)) { // Set b0 extraction as fixed image - b0Extraction->SetInput(fixedDwi->GetVectorImage()); - b0Extraction->SetDirections(fixedDwi->GetDirections()); + b0Extraction->SetInput( itkFixedDwiPointer ); + b0Extraction->SetDirections(fixedGradientDirections); b0Extraction->Update(); mitk::Image::Pointer tmp = mitk::Image::New(); tmp->InitializeByItk(b0Extraction->GetOutput()); tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer()); fixedImage = tmp; } - if (movingDwi != NULL) + if (movingGradientDirections || ( movingGradientDirections->Size() == 0)) { // Set b0 extraction as moving image - b0Extraction->SetInput(movingDwi->GetVectorImage()); - b0Extraction->SetDirections(movingDwi->GetDirections()); + b0Extraction->SetInput( itkMovingDwiPointer ); + b0Extraction->SetDirections(movingGradientDirections); b0Extraction->Update(); mitk::Image::Pointer tmp = mitk::Image::New(); tmp->InitializeByItk(b0Extraction->GetOutput()); tmp->SetVolume(b0Extraction->GetOutput()->GetBufferPointer()); movingImage = tmp; } // align the offsets of the two images. this is done to avoid non-overlapping initialization Point3D origin = fixedImage->GetGeometry()->GetOrigin(); Point3D originMoving = movingImage->GetGeometry()->GetOrigin(); mitk::Image::Pointer tmpImage = movingImage->Clone(); if (useSameOrigin) { offset[0] = originMoving[0]-origin[0]; offset[1] = originMoving[1]-origin[1]; offset[2] = originMoving[2]-origin[2]; tmpImage->GetGeometry()->SetOrigin(origin); } // Start registration mitk::PyramidImageRegistrationMethod::Pointer registrationMethod = mitk::PyramidImageRegistrationMethod::New(); registrationMethod->SetFixedImage( fixedImage ); if (mask != NULL) { registrationMethod->SetFixedImageMask(mask); registrationMethod->SetUseFixedImageMask(true); } else { registrationMethod->SetUseFixedImageMask(false); } registrationMethod->SetTransformToRigid(); registrationMethod->SetCrossModalityOn(); registrationMethod->SetMovingImage(tmpImage); registrationMethod->Update(); registrationMethod->GetParameters(transformation); // first three: euler angles, last three translation } diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h index deead178e9..7b7802c3cc 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkAdcImageFilter.h @@ -1,80 +1,80 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkAdcImageFilter_h_ #define __itkAdcImageFilter_h_ #include "itkImageToImageFilter.h" #include "itkVectorImage.h" -#include +#include namespace itk{ /** \class AdcImageFilter */ template< class TInPixelType, class TOutPixelType > class AdcImageFilter : public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > { public: typedef AdcImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > Superclass; - typedef mitk::DiffusionImage< short >::GradientDirectionType GradientDirectionType; - typedef mitk::DiffusionImage< short >::GradientDirectionContainerType::Pointer GradientContainerType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GradientContainerType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(AdcImageFilter, ImageToImageFilter) typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; itkSetMacro( B_value, double ) itkSetMacro( GradientDirections, GradientContainerType ) protected: AdcImageFilter(); ~AdcImageFilter() {} void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); double m_B_value; GradientContainerType m_GradientDirections; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkAdcImageFilter.txx" #endif #endif //__itkAdcImageFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h index cdb43eec2e..c268bb5e40 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDwiNormilzationFilter.h @@ -1,100 +1,100 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkDwiNormilzationFilter_h_ #define __itkDwiNormilzationFilter_h_ #include "itkImageToImageFilter.h" #include "itkVectorImage.h" -#include #include #include #include +#include namespace itk{ /** \class DwiNormilzationFilter * \brief Normalizes the data vectors either using the specified reference value or the voxelwise baseline value. */ template< class TInPixelType > class DwiNormilzationFilter : public ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > { public: typedef DwiNormilzationFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(DwiNormilzationFilter, ImageToImageFilter) typedef itk::Image< double, 3 > DoubleImageType; typedef itk::Image< unsigned char, 3 > UcharImageType; typedef itk::Image< unsigned short, 3 > BinImageType; typedef itk::Image< TInPixelType, 3 > TInPixelImageType; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - typedef mitk::DiffusionImage< short >::GradientDirectionType GradientDirectionType; - typedef mitk::DiffusionImage< short >::GradientDirectionContainerType::Pointer GradientContainerType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GradientContainerType; typedef itk::LabelStatisticsImageFilter< TInPixelImageType,BinImageType > StatisticsFilterType; typedef itk::Statistics::Histogram< typename TInPixelImageType::PixelType > HistogramType; typedef typename HistogramType::MeasurementType MeasurementType; typedef itk::ShiftScaleImageFilter ShiftScaleImageFilterType; itkSetMacro( GradientDirections, GradientContainerType ) itkSetMacro( ScalingFactor, TInPixelType ) itkSetMacro( UseGlobalReference, bool ) itkSetMacro( MaskImage, UcharImageType::Pointer ) itkSetMacro( Reference, double ) protected: DwiNormilzationFilter(); ~DwiNormilzationFilter() {} void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); UcharImageType::Pointer m_MaskImage; GradientContainerType m_GradientDirections; int m_B0Index; TInPixelType m_ScalingFactor; bool m_UseGlobalReference; double m_Reference; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDwiNormilzationFilter.txx" #endif #endif //__itkDwiNormilzationFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h index dd75de3f14..9c80bb27e4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h @@ -1,123 +1,126 @@ /*=================================================================== 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/itkElectrostaticRepulsionDiffusionGradientReductionFilter.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_ElectrostaticRepulsionDiffusionGradientReductionFilter_h_ #define _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_h_ #include #include #include namespace itk { /** * \brief Select subset of the input vectors equally distributed over the sphere using an iterative electrostatic repulsion strategy. */ template class ElectrostaticRepulsionDiffusionGradientReductionFilter : public ImageToImageFilter, itk::VectorImage > { public: typedef ElectrostaticRepulsionDiffusionGradientReductionFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< itk::VectorImage, itk::VectorImage > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(ElectrostaticRepulsionDiffusionGradientReductionFilter, ImageToImageFilter) typedef TInputScalarType InputScalarType; typedef itk::VectorImage InputImageType; typedef typename InputImageType::PixelType InputPixelType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef std::vector IndicesVector; typedef std::map BValueMap; itkGetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkSetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkGetMacro(GradientDirections, GradientDirectionContainerType::Pointer) itkSetMacro(GradientDirections, GradientDirectionContainerType::Pointer) IndicesVector GetUsedGradientIndices(){return m_UsedGradientIndices;} void SetOriginalBValueMap(BValueMap inp){m_OriginalBValueMap = inp;} void SetShellSelectionBValueMap(BValueMap inp){m_InputBValueMap = inp;} void SetNumGradientDirections(std::vector numDirs){m_NumGradientDirections = numDirs;} + + void UpdateOutputInformation(); + protected: ElectrostaticRepulsionDiffusionGradientReductionFilter(); ~ElectrostaticRepulsionDiffusionGradientReductionFilter() {} void GenerateData(); double Costs(); ///< calculates electrostatic energy of current direction set GradientDirectionContainerType::Pointer m_GradientDirections; ///< container for the subsampled output gradient directions GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions IndicesVector m_UsedGradientIndices; IndicesVector m_UnusedGradientIndices; IndicesVector m_BaselineImageIndices; BValueMap m_OriginalBValueMap; BValueMap m_InputBValueMap; std::vector m_NumGradientDirections; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx index 3b601cdb0a..7253e35b42 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx @@ -1,247 +1,261 @@ /*=================================================================== 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/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_txx_ #define _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_txx_ #endif #define _USE_MATH_DEFINES #include "itkElectrostaticRepulsionDiffusionGradientReductionFilter.h" #include #include #include #include namespace itk { template ElectrostaticRepulsionDiffusionGradientReductionFilter ::ElectrostaticRepulsionDiffusionGradientReductionFilter() { this->SetNumberOfRequiredInputs( 1 ); } template double ElectrostaticRepulsionDiffusionGradientReductionFilter ::Costs() { double costs = 2*M_PI; for (IndicesVector::iterator it = m_UsedGradientIndices.begin(); it!=m_UsedGradientIndices.end(); ++it) { for (IndicesVector::iterator it2 = m_UsedGradientIndices.begin(); it2!=m_UsedGradientIndices.end(); ++it2) if (it != it2) { vnl_vector_fixed v1 = m_OriginalGradientDirections->at(*it); vnl_vector_fixed v2 = m_OriginalGradientDirections->at(*it2); v1.normalize(); v2.normalize(); double temp = dot_product(v1,v2); if (temp>1) temp = 1; else if (temp<-1) temp = -1; double angle = acos(temp); if (angle1) temp = 1; else if (temp<-1) temp = -1; angle = acos(temp); if (angle void ElectrostaticRepulsionDiffusionGradientReductionFilter ::GenerateData() { unsigned int randSeed = time(NULL); if(m_InputBValueMap.empty() || m_NumGradientDirections.size()!=m_InputBValueMap.size()) return; BValueMap manipulatedMap = m_InputBValueMap; int shellCounter = 0; for(BValueMap::iterator it = m_InputBValueMap.begin(); it != m_InputBValueMap.end(); it++ ) { srand(randSeed); // initialize index vectors m_UsedGradientIndices.clear(); m_UnusedGradientIndices.clear(); if ( it->second.size() <= m_NumGradientDirections[shellCounter] ) { itkWarningMacro( << "current directions: " << it->second.size() << " wanted directions: " << m_NumGradientDirections[shellCounter]); m_NumGradientDirections[shellCounter] = it->second.size(); shellCounter++; continue; } MITK_INFO << "Shell number: " << shellCounter; unsigned int c=0; for (unsigned int i=0; isecond.size(); i++) { if (csecond.at(i)); else m_UnusedGradientIndices.push_back(it->second.at(i)); c++; } double minAngle = Costs(); double newMinAngle = 0; MITK_INFO << "minimum angle: " << 180*minAngle/M_PI; int stagnationCount = 0; int rejectionCount = 0; int maxRejections = m_NumGradientDirections[shellCounter] * 1000; if (maxRejections<10000) maxRejections = 10000; int iUsed = 0; if (m_UsedGradientIndices.size()>0) while ( stagnationCount<1000 && rejectionCount minAngle) // accept or reject proposal { MITK_INFO << "minimum angle: " << 180*newMinAngle/M_PI; if ( (newMinAngle-minAngle)<0.01 ) stagnationCount++; else stagnationCount = 0; minAngle = newMinAngle; rejectionCount = 0; } else { rejectionCount++; m_UsedGradientIndices.at(iUsed) = vUsed; m_UnusedGradientIndices.at(iUnUsed) = vUnUsed; } iUsed++; iUsed = iUsed % m_UsedGradientIndices.size(); } manipulatedMap[it->first] = m_UsedGradientIndices; shellCounter++; } int vecLength = 0 ; for(BValueMap::iterator it = manipulatedMap.begin(); it != manipulatedMap.end(); it++) vecLength += it->second.size(); // initialize output image 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( vecLength ); // Set the vector length outImage->Allocate(); itk::ImageRegionIterator< OutputImageType > newIt(outImage, outImage->GetLargestPossibleRegion()); newIt.GoToBegin(); typename InputImageType::Pointer inImage = const_cast(this->GetInput(0)); itk::ImageRegionIterator< InputImageType > oldIt(inImage, inImage->GetLargestPossibleRegion()); oldIt.GoToBegin(); // initial new value of voxel OutputPixelType newVec; newVec.SetSize( vecLength ); newVec.AllocateElements( vecLength ); // generate new pixel values while(!newIt.IsAtEnd()) { // init new vector with zeros newVec.Fill(0.0); InputPixelType oldVec = oldIt.Get(); int index = 0; for(BValueMap::iterator it=manipulatedMap.begin(); it!=manipulatedMap.end(); it++) for(unsigned int j=0; jsecond.size(); j++) { newVec[index] = oldVec[it->second.at(j)]; index++; } newIt.Set(newVec); ++newIt; ++oldIt; } // set new gradient directions m_GradientDirections = GradientDirectionContainerType::New(); int index = 0; for(BValueMap::iterator it = manipulatedMap.begin(); it != manipulatedMap.end(); it++) for(unsigned int j = 0; j < it->second.size(); j++) { m_GradientDirections->InsertElement(index, m_OriginalGradientDirections->at(it->second.at(j))); index++; } this->SetNumberOfRequiredOutputs (1); this->SetNthOutput (0, outImage); MITK_INFO << "...done"; } +template +void +ElectrostaticRepulsionDiffusionGradientReductionFilter +::UpdateOutputInformation() +{ + Superclass::UpdateOutputInformation(); + int vecLength = 0 ; + for(unsigned int index = 0; index < m_NumGradientDirections.size(); index++) + { + vecLength += m_NumGradientDirections[index]; + } + + this->GetOutput()->SetVectorLength( vecLength ); +} } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h index 5db210b559..54fcc4e434 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkExtractDwiChannelFilter.h @@ -1,78 +1,77 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkExtractDwiChannelFilter_h_ #define __itkExtractDwiChannelFilter_h_ #include "itkImageToImageFilter.h" #include "itkVectorImage.h" -#include #include -#include +#include namespace itk{ /** \class ExtractDwiChannelFilter - * \brief Remove spcified channels from diffusion-weighted image. + * \brief Remove specified channels from diffusion-weighted image. */ template< class TInPixelType > class ExtractDwiChannelFilter : public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TInPixelType, 3 > > { public: typedef ExtractDwiChannelFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TInPixelType, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(ExtractDwiChannelFilter, ImageToImageFilter) typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; itkSetMacro( ChannelIndex, unsigned int ) protected: ExtractDwiChannelFilter(); ~ExtractDwiChannelFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType ); unsigned int m_ChannelIndex; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkExtractDwiChannelFilter.txx" #endif #endif //__itkExtractDwiChannelFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h index 878443938b..34cf4ef219 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h @@ -1,147 +1,145 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkNonLocalMeansDenoisingFilter_h_ #define __itkNonLocalMeansDenoisingFilter_h_ #include "itkImageToImageFilter.h" #include "itkVectorImage.h" -#include - namespace itk{ /** @class NonLocalMeansDenoisingFilter * @brief This class denoises a vectorimage according to the non-local means procedure. * * This Filter needs as an input a diffusion weigthed image, which will be denoised unsing the non-local means principle. * An input mask is optional to denoise only inside the mask range. All other voxels will be set to 0. */ template< class TPixelType > class NonLocalMeansDenoisingFilter : public ImageToImageFilter< VectorImage < TPixelType, 3 >, VectorImage < TPixelType, 3 > > { public: /** Typedefs */ typedef NonLocalMeansDenoisingFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage < TPixelType, 3 >, VectorImage < TPixelType, 3 > > Superclass; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef Image MaskImageType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(NonLocalMeansDenoisingFilter, ImageToImageFilter) /** * @brief Set flag to use joint information */ itkSetMacro(UseJointInformation, bool) /** * @brief Set the searchradius * * The searchradius generates a neighborhood of size (2 * searchradius + 1)³. * Default is 4. */ itkSetMacro(SearchRadius, int) /** * @brief Set the comparisonradius * * The comparisonradius generates neighborhoods of size (2 * comparisonradius +1)³. * Default is 1. */ itkSetMacro(ComparisonRadius, int) /** * @brief Set the variance of the noise * * The variance of the noise needs to be estimated to use this filter properly. * Default is 1. */ itkSetMacro(Variance, double) /** * @brief Set flag to use a rician adaption * * If this flag is true the filter uses a method which is optimized for Rician distributed noise. */ itkSetMacro(UseRicianAdaption, bool) /** * @brief Get the amount of calculated Voxels * * @return the number of calculated Voxels until yet, useful for the use of a progressbars. */ itkGetMacro(CurrentVoxelCount, unsigned int) /** @brief Set the input image. **/ void SetInputImage(const InputImageType* image); /** * @brief Set a denoising mask * * optional * * Set a mask to denoise only the masked area, all voxel outside this area will be set to 0. */ void SetInputMask(MaskImageType* mask); protected: NonLocalMeansDenoisingFilter(); ~NonLocalMeansDenoisingFilter() {} /** * @brief Calculations which need to be done before the denoising starts * * This method is called before the denoising starts. It calculates the ROI if a mask is used * and sets the number of processed voxels to zero. */ void BeforeThreadedGenerateData(); /** * @brief Denoising procedure * * This method calculates the denoised voxelvalue for each voxel in the image in multiple threads. * If a mask is used, voxels outside the masked area will be set to 0. * * @param outputRegionForThread Region to denoise for each thread. */ void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); private: int m_SearchRadius; ///< Radius of the searchblock. int m_ComparisonRadius; ///< Radius of the comparisonblock. bool m_UseJointInformation; ///< Flag to use joint information. bool m_UseRicianAdaption; ///< Flag to use rician adaption. unsigned int m_CurrentVoxelCount; ///< Amount of processed voxels. double m_Variance; ///< Estimated noise variance. typename MaskImageType::Pointer m_Mask; ///< Pointer to the mask image. }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkNonLocalMeansDenoisingFilter.txx" #endif #endif //__itkNonLocalMeansDenoisingFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h index 7e18435d36..168a92b5f8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkRemoveDwiChannelFilter.h @@ -1,85 +1,84 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkRemoveDwiChannelFilter_h_ #define __itkRemoveDwiChannelFilter_h_ #include "itkImageToImageFilter.h" #include "itkVectorImage.h" -#include #include -#include +#include namespace itk{ /** \class RemoveDwiChannelFilter * \brief Remove spcified channels from diffusion-weighted image. */ template< class TInPixelType > class RemoveDwiChannelFilter : public ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > { public: typedef RemoveDwiChannelFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, VectorImage< TInPixelType, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(RemoveDwiChannelFilter, ImageToImageFilter) typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - typedef typename mitk::DiffusionImage< TInPixelType >::GradientDirectionType DirectionType; - typedef typename mitk::DiffusionImage< TInPixelType >::GradientDirectionContainerType DirectionContainerType; + typedef typename mitk::DiffusionPropertyHelper::GradientDirectionType DirectionType; + typedef typename mitk::DiffusionPropertyHelper::GradientDirectionsContainerType DirectionContainerType; void SetChannelIndices( std::vector< unsigned int > indices ){ m_ChannelIndices = indices; } void SetDirections( typename DirectionContainerType::Pointer directions ){ m_Directions = directions; } typename DirectionContainerType::Pointer GetNewDirections(){ return m_NewDirections; } protected: RemoveDwiChannelFilter(); ~RemoveDwiChannelFilter() {} void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType id ); std::vector< unsigned int > m_ChannelIndices; typename DirectionContainerType::Pointer m_Directions; typename DirectionContainerType::Pointer m_NewDirections; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkRemoveDwiChannelFilter.txx" #endif #endif //__itkRemoveDwiChannelFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h index 2bbaf13288..925d1ff77a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.h @@ -1,138 +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. ===================================================================*/ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkResampleDwiImageFilter.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_ResampleDwiImageFilter_h_ #define _itk_ResampleDwiImageFilter_h_ #include #include #include namespace itk { /** * \brief Resample DWI channel by channel. */ template class ResampleDwiImageFilter : public ImageToImageFilter, itk::VectorImage > { public: enum Interpolation { Interpolate_NearestNeighbour, Interpolate_Linear, Interpolate_BSpline, Interpolate_WindowedSinc }; typedef ResampleDwiImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< itk::VectorImage, itk::VectorImage > Superclass; typedef itk::Vector< double, 3 > DoubleVectorType; typedef itk::VectorImage DwiImageType; typedef itk::Image DwiChannelType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(ResampleDwiImageFilter, ImageToImageFilter) itkSetMacro( Interpolation, Interpolation ) void SetSamplingFactor(DoubleVectorType sampling) { m_NewSpacing = this->GetInput()->GetSpacing(); m_NewSpacing[0] /= sampling[0]; m_NewSpacing[1] /= sampling[1]; m_NewSpacing[2] /= sampling[2]; m_NewImageRegion = this->GetInput()->GetLargestPossibleRegion(); m_NewImageRegion.SetSize(0, m_NewImageRegion.GetSize(0)*sampling[0]); m_NewImageRegion.SetSize(1, m_NewImageRegion.GetSize(1)*sampling[1]); m_NewImageRegion.SetSize(2, m_NewImageRegion.GetSize(2)*sampling[2]); } void SetNewSpacing(DoubleVectorType spacing) { DoubleVectorType oldSpacing = this->GetInput()->GetSpacing(); DoubleVectorType sampling; sampling[0] = oldSpacing[0]/spacing[0]; sampling[1] = oldSpacing[1]/spacing[1]; sampling[2] = oldSpacing[2]/spacing[2]; m_NewSpacing = spacing; m_NewImageRegion = this->GetInput()->GetLargestPossibleRegion(); m_NewImageRegion.SetSize(0, m_NewImageRegion.GetSize(0)*sampling[0]); m_NewImageRegion.SetSize(1, m_NewImageRegion.GetSize(1)*sampling[1]); m_NewImageRegion.SetSize(2, m_NewImageRegion.GetSize(2)*sampling[2]); } void SetNewImageSize(ImageRegion<3> region) { ImageRegion<3> oldRegion = this->GetInput()->GetLargestPossibleRegion(); DoubleVectorType sampling; sampling[0] = (double)region.GetSize(0)/oldRegion.GetSize(0); sampling[1] = (double)region.GetSize(1)/oldRegion.GetSize(1); sampling[2] = (double)region.GetSize(2)/oldRegion.GetSize(2); m_NewImageRegion = region; m_NewSpacing = this->GetInput()->GetSpacing(); m_NewSpacing[0] /= sampling[0]; m_NewSpacing[1] /= sampling[1]; m_NewSpacing[2] /= sampling[2]; } + virtual void UpdateOutputInformation(); + + virtual void PropagateRequestedRegion(){} + + virtual void PropagateRequestedRegion(itk::DataObject *output){} + + virtual void VerifyInputInformation(){} + protected: ResampleDwiImageFilter(); ~ResampleDwiImageFilter(){} void GenerateData(); DoubleVectorType m_NewSpacing; ImageRegion<3> m_NewImageRegion; Interpolation m_Interpolation; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkResampleDwiImageFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx index 516c121104..b317f88f4d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResampleDwiImageFilter.txx @@ -1,166 +1,176 @@ /*=================================================================== 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/itkResampleDwiImageFilter.txx $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _itk_ResampleDwiImageFilter_txx_ #define _itk_ResampleDwiImageFilter_txx_ #endif #define _USE_MATH_DEFINES #include "itkResampleDwiImageFilter.h" #include #include #include #include #include #include #include namespace itk { template ResampleDwiImageFilter ::ResampleDwiImageFilter() : m_Interpolation(Interpolate_Linear) { this->SetNumberOfRequiredInputs( 1 ); } template void ResampleDwiImageFilter ::GenerateData() { // // initialize output image // itk::Vector< double, 3 > spacing = this->GetInput()->GetSpacing(); // spacing[0] /= m_SamplingFactor[0]; // spacing[1] /= m_SamplingFactor[1]; // spacing[2] /= m_SamplingFactor[2]; // ImageRegion<3> region = this->GetInput()->GetLargestPossibleRegion(); // region.SetSize(0, region.GetSize(0)*m_SamplingFactor[0]); // region.SetSize(1, region.GetSize(1)*m_SamplingFactor[1]); // region.SetSize(2, region.GetSize(2)*m_SamplingFactor[2]); itk::Point origin = this->GetInput()->GetOrigin(); origin[0] -= this->GetInput()->GetSpacing()[0]/2; origin[1] -= this->GetInput()->GetSpacing()[1]/2; origin[2] -= this->GetInput()->GetSpacing()[2]/2; origin[0] += m_NewSpacing[0]/2; origin[1] += m_NewSpacing[1]/2; origin[2] += m_NewSpacing[2]/2; typename DwiImageType::Pointer outImage = DwiImageType::New(); outImage->SetSpacing( m_NewSpacing ); outImage->SetOrigin( origin ); outImage->SetDirection( this->GetInput()->GetDirection() ); outImage->SetLargestPossibleRegion( m_NewImageRegion ); outImage->SetBufferedRegion( m_NewImageRegion ); outImage->SetRequestedRegion( m_NewImageRegion ); outImage->SetVectorLength( this->GetInput()->GetVectorLength() ); outImage->Allocate(); typename itk::ResampleImageFilter::Pointer resampler = itk::ResampleImageFilter::New(); resampler->SetOutputParametersFromImage(outImage); switch (m_Interpolation) { case Interpolate_NearestNeighbour: { typename itk::NearestNeighborInterpolateImageFunction::Pointer interp = itk::NearestNeighborInterpolateImageFunction::New(); resampler->SetInterpolator(interp); break; } case Interpolate_Linear: { typename itk::LinearInterpolateImageFunction::Pointer interp = itk::LinearInterpolateImageFunction::New(); resampler->SetInterpolator(interp); break; } case Interpolate_BSpline: { typename itk::BSplineInterpolateImageFunction::Pointer interp = itk::BSplineInterpolateImageFunction::New(); resampler->SetInterpolator(interp); break; } case Interpolate_WindowedSinc: { typename itk::WindowedSincInterpolateImageFunction::Pointer interp = itk::WindowedSincInterpolateImageFunction::New(); resampler->SetInterpolator(interp); break; } default: { typename itk::LinearInterpolateImageFunction::Pointer interp = itk::LinearInterpolateImageFunction::New(); resampler->SetInterpolator(interp); } } for (unsigned int i=0; iGetInput()->GetVectorLength(); i++) { typename DwiChannelType::Pointer channel = DwiChannelType::New(); channel->SetSpacing( this->GetInput()->GetSpacing() ); channel->SetOrigin( this->GetInput()->GetOrigin() ); channel->SetDirection( this->GetInput()->GetDirection() ); channel->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion() ); channel->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() ); channel->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() ); channel->Allocate(); ImageRegionIterator it(channel, channel->GetLargestPossibleRegion()); while(!it.IsAtEnd()) { typename DwiImageType::PixelType pix = this->GetInput()->GetPixel(it.GetIndex()); it.Set(pix[i]); ++it; } resampler->SetInput(channel); resampler->Update(); channel = resampler->GetOutput(); ImageRegionIterator it2(outImage, outImage->GetLargestPossibleRegion()); while(!it2.IsAtEnd()) { typename DwiImageType::PixelType pix = it2.Get(); pix[i] = channel->GetPixel(it2.GetIndex()); it2.Set(pix); ++it2; } } this->SetNthOutput(0, outImage); } +template +void +ResampleDwiImageFilter +::UpdateOutputInformation() +{ + // Calls to superclass updateoutputinformation + //Superclass::UpdateOutputInformation(); + this->GetOutput()->SetSpacing(m_NewSpacing); + this->GetOutput()->SetLargestPossibleRegion(m_NewImageRegion); +} } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx index 49715a1226..fb1c06a99a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkResidualImageFilter.txx @@ -1,285 +1,285 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.txx $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _itk_ResidualImageFilter_txx_ #define _itk_ResidualImageFilter_txx_ #endif #include "itkResidualImageFilter.h" #include #include namespace itk { template void ResidualImageFilter ::GenerateData() { typename InputImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize(); typename InputImageType::SizeType size2 = m_SecondDiffusionImage->GetLargestPossibleRegion().GetSize(); if(size != size2) { MITK_ERROR << "Sizes do not match"; return; } // Initialize output image typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); outputImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing outputImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin outputImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction outputImage->SetRegions( this->GetInput()->GetLargestPossibleRegion() ); outputImage->Allocate(); outputImage->FillBuffer(0.0); std::vector< std::vector > residuals; // per slice, per volume std::vector< std::vector > > residualsPerSlice; // Detrmine number of B0 images int numberB0=0; for(unsigned int i=0; iSize(); i++) { GradientDirectionType grad = m_Gradients->ElementAt(i); if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001) { numberB0++; } } residuals.resize(this->GetInput()->GetVectorLength()-numberB0); // Calculate the standard residual image and for each volume put all residuals in a vector for(unsigned int z=0; z > sliceResiduals; // residuals per volume for this slice sliceResiduals.resize(this->GetInput()->GetVectorLength()-numberB0); for(unsigned int y=0; y ix; ix[0] = x; ix[1] = y; ix[2] = z; typename InputImageType::PixelType p1 = this->GetInput()->GetPixel(ix); typename InputImageType::PixelType p2 = m_SecondDiffusionImage->GetPixel(ix); if(p1.GetSize() != p2.GetSize()) { MITK_ERROR << "Vector sizes do not match"; return; } if(p1.GetElement(m_B0Index) <= m_B0Threshold) { continue; } double res = 0; int shift = 0; // correction for the skipped B0 images for(unsigned int i = 0; iElementAt(i); if(!(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001)) { double val1 = (double)p1.GetElement(i); double val2 = (double)p2.GetElement(i); res += abs(val1-val2); residuals[i-shift].push_back(val1-val2); sliceResiduals[i-shift].push_back(val1-val2); } else { shift++; } } res = res/p1.GetSize(); outputImage->SetPixel(ix, res); } } residualsPerSlice.push_back(sliceResiduals); } // for each dw volume: sort the the measured residuals (for each voxel) to enable determining Q1 and Q3; calculate means // determine percentage of errors as described in QUALITY ASSESSMENT THROUGH ANALYSIS OF RESIDUALS OF DIFFUSION IMAGE FITTING // Leemans et al 2008 double q1,q3, median; std::vector< std::vector >::iterator it = residuals.begin(); while(it != residuals.end()) { std::vector res = *it; // sort std::sort(res.begin(), res.end()); q1 = res[0.25*res.size()]; m_Q1.push_back(q1); q3 = res[0.75*res.size()]; m_Q3.push_back(q3); median = res[0.5*res.size()]; double iqr = q3-q1; double outlierThreshold = median + 1.5*iqr; double numberOfOutliers = 0.0; std::vector::iterator resIt = res.begin(); double mean = 0; while(resIt != res.end()) { double f = *resIt; if(f>outlierThreshold) { numberOfOutliers++; } mean += f; ++resIt; } double percOfOutliers = 100 * numberOfOutliers / res.size(); m_PercentagesOfOutliers.push_back(percOfOutliers); mean /= res.size(); m_Means.push_back(mean); ++it; } // Calculate for each slice the number of outliers per volume(dw volume) std::vector< std::vector > >::iterator sliceIt = residualsPerSlice.begin(); while(sliceIt != residualsPerSlice.end()) { std::vector< std::vector > currentSlice = *sliceIt; std::vector percentages; std::vector< std::vector >::iterator volIt = currentSlice.begin(); while(volIt != currentSlice.end()) { std::vector currentVolume = *volIt; //sort std::sort(currentVolume.begin(), currentVolume.end()); q1 = currentVolume[0.25*currentVolume.size()]; q3 = currentVolume[0.75*currentVolume.size()]; median = currentVolume[0.5*currentVolume.size()]; double iqr = q3-q1; double outlierThreshold = median + 1.5*iqr; double numberOfOutliers = 0.0; std::vector::iterator resIt = currentVolume.begin(); - double mean; + double mean(0.0); while(resIt != currentVolume.end()) { double f = *resIt; if(f>outlierThreshold) { numberOfOutliers++; } mean += f; ++resIt; } double percOfOutliers = 100 * numberOfOutliers / currentVolume.size(); percentages.push_back(percOfOutliers); ++volIt; } m_OutliersPerSlice.push_back(percentages); ++sliceIt; } } } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx index 34cb1f08cb..ad119b1a73 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorDerivedMeasurementsFilter.txx @@ -1,145 +1,145 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkTensorDerivedMeasurementsFilter_txx #define __itkTensorDerivedMeasurementsFilter_txx - +#include namespace itk { template TensorDerivedMeasurementsFilter::TensorDerivedMeasurementsFilter() : m_Measure(AD) { } template void TensorDerivedMeasurementsFilter::GenerateData() { typename TensorImageType::Pointer tensorImage = static_cast< TensorImageType * >( this->ProcessObject::GetInput(0) ); typedef ImageRegionConstIterator< TensorImageType > TensorImageIteratorType; typedef ImageRegionIterator< OutputImageType > OutputImageIteratorType; typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); typename TensorImageType::RegionType region = tensorImage->GetLargestPossibleRegion(); outputImage->SetSpacing( tensorImage->GetSpacing() ); // Set the image spacing outputImage->SetOrigin( tensorImage->GetOrigin() ); // Set the image origin outputImage->SetDirection( tensorImage->GetDirection() ); // Set the image direction outputImage->SetRegions( tensorImage->GetLargestPossibleRegion()); outputImage->Allocate(); TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetLargestPossibleRegion()); OutputImageIteratorType outputIt(outputImage, outputImage->GetLargestPossibleRegion()); tensorIt.GoToBegin(); outputIt.GoToBegin(); while(!tensorIt.IsAtEnd() && !outputIt.IsAtEnd()){ TensorType tensor = tensorIt.Get(); switch(m_Measure) { case FA: { TPixel diffusionIndex = tensor.GetFractionalAnisotropy(); outputIt.Set(diffusionIndex); break; } case RA: { TPixel diffusionIndex = tensor.GetRelativeAnisotropy(); outputIt.Set(diffusionIndex); break; } case AD: { // eigenvalues are sorted in ascending order by default because the // itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation typename TensorType::EigenValuesArrayType evs; tensor.ComputeEigenValues(evs); outputIt.Set(evs[2]); break; } case RD: { // eigenvalues are sorted in ascending order by default because the // itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation typename TensorType::EigenValuesArrayType evs; tensor.ComputeEigenValues(evs); outputIt.Set((evs[0]+evs[1])/2.0); break; } case CA: { // eigenvalues are sorted in ascending order by default because the // itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation typename TensorType::EigenValuesArrayType evs; tensor.ComputeEigenValues(evs); if (evs[2] == 0) { outputIt.Set(0); break; } outputIt.Set(1.0-(evs[0]+evs[1])/(2.0*evs[2])); break; } case L2: { // eigenvalues are sorted in ascending order by default because the // itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation typename TensorType::EigenValuesArrayType evs; tensor.ComputeEigenValues(evs); outputIt.Set(evs[1]); break; } case L3: { // eigenvalues are sorted in ascending order by default because the // itk::SymmetricEigenAnalysis defaults are not touched in the tensor implementation typename TensorType::EigenValuesArrayType evs; tensor.ComputeEigenValues(evs); outputIt.Set(evs[0]); break; } case MD: { typename TensorType::EigenValuesArrayType evs; tensor.ComputeEigenValues(evs); outputIt.Set((evs[0]+evs[1]+evs[2])/3.0); break; } } ++tensorIt; ++outputIt; } } } #endif // __itkTensorDerivedMeasurements_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h index cbed2d792b..b588be572a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.h @@ -1,168 +1,171 @@ /*=================================================================== 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 namespace itk { template class TensorImageToDiffusionImageFilter : public ImageToImageFilter,3>, itk::VectorImage > { public: typedef TInputScalarType InputScalarType; typedef itk::DiffusionTensor3D InputPixelType; typedef itk::Image InputImageType; typedef typename InputImageType::RegionType InputImageRegionType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef typename OutputImageType::RegionType OutputImageRegionType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef typename BaselineImageType::RegionType BaselineImageRegionType; typedef itk::Image< short, 3> MaskImageType; typedef MaskImageType::RegionType MaskImageRegionType; typedef TensorImageToDiffusionImageFilter Self; typedef ImageToImageFilter Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; itkTypeMacro (TensorImageToDiffusionImageFilter, ImageToImageFilter); itkStaticConstMacro (ImageDimension, unsigned int, OutputImageType::ImageDimension); itkFactorylessNewMacro(Self) itkCloneMacro(Self) typedef vnl_vector_fixed GradientType; typedef VectorContainer GradientListType; typedef GradientListType::Pointer GradientListPointerType; /** Manually Set/Get a list of gradients */ void SetGradientList(const GradientListPointerType list) { m_GradientList = list; this->Modified(); } GradientListPointerType GetGradientList(void) const {return m_GradientList;} void SetBValue( const double& bval) { m_BValue = bval; } void SetMaskImage( MaskImageType::Pointer maskimage ) { m_MaskImage = maskimage; } /** * @brief Set an external baseline image for signal generation (optional) * * An option to enforce a specific baseline image. If none provided (default) the filter uses * the itk::TensorToL2NormImageFilter to generate the modelled baseline image. */ void SetExternalBaselineImage( typename BaselineImageType::Pointer bimage) { m_BaselineImage = bimage; } itkSetMacro(Min, OutputScalarType); itkSetMacro(Max, OutputScalarType); protected: TensorImageToDiffusionImageFilter() { m_BValue = 1.0; m_BaselineImage = 0; m_Min = 0.0; m_Max = 10000.0; } virtual ~TensorImageToDiffusionImageFilter(){} void PrintSelf (std::ostream& os, Indent indent) const { Superclass::PrintSelf (os, indent); } virtual void BeforeThreadedGenerateData( void ); virtual void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); //void GenerateData(); - GradientListPointerType m_GradientList; + virtual void UpdateOutputInformation(); + + + private: + + TensorImageToDiffusionImageFilter (const Self&); + void operator=(const Self&); + + GradientListType::Pointer m_GradientList; double m_BValue; typename BaselineImageType::Pointer m_BaselineImage; OutputScalarType m_Min; OutputScalarType m_Max; MaskImageType::Pointer m_MaskImage; - private: - - TensorImageToDiffusionImageFilter (const Self&); - void operator=(const Self&); - }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkTensorImageToDiffusionImageFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx index e812dfb51a..67a9b355d2 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkTensorImageToDiffusionImageFilter.txx @@ -1,235 +1,244 @@ /*=================================================================== 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 ::ThreadedGenerateData (const OutputImageRegionType &outputRegionForThread, ThreadIdType threadId ) { typedef ImageRegionIterator IteratorOutputType; typedef ImageRegionConstIterator IteratorInputType; typedef ImageRegionConstIterator IteratorBaselineType; unsigned long numPixels = outputRegionForThread.GetNumberOfPixels(); unsigned long step = numPixels/100; unsigned long progress = 0; IteratorOutputType itOut (this->GetOutput(), outputRegionForThread); IteratorInputType itIn (this->GetInput(), outputRegionForThread); IteratorBaselineType itB0 (m_BaselineImage, outputRegionForThread); typedef ImageRegionConstIterator< MaskImageType > IteratorMaskImageType; IteratorMaskImageType itMask; if( m_MaskImage.IsNotNull() ) { itMask = IteratorMaskImageType( m_MaskImage, outputRegionForThread); itMask.GoToBegin(); } if( threadId==0 ) { this->UpdateProgress (0.0); } while(!itIn.IsAtEnd()) { if( this->GetAbortGenerateData() ) { throw itk::ProcessAborted(__FILE__,__LINE__); } InputPixelType T = itIn.Get(); BaselinePixelType b0 = itB0.Get(); OutputPixelType out; out.SetSize(m_GradientList->Size()); out.Fill(0); short maskvalue = 1; if( m_MaskImage.IsNotNull() ) { maskvalue = itMask.Get(); ++itMask; } std::vector b0_indices; if( b0 > 0) { for( unsigned int i=0; iSize(); i++) { GradientType g = m_GradientList->at(i); // normalize vector so the following computations work const double twonorm = g.two_norm(); if( twonorm < vnl_math::eps ) { b0_indices.push_back(i); continue; } GradientType gn = g.normalize(); InputPixelType S; S[0] = gn[0]*gn[0]; S[1] = gn[1]*gn[0]; S[2] = gn[2]*gn[0]; S[3] = gn[1]*gn[1]; S[4] = gn[2]*gn[1]; S[5] = gn[2]*gn[2]; const double res = T[0]*S[0] + 2 * T[1]*S[1] + T[3]*S[3] + 2 * T[2]*S[2] + 2 * T[4]*S[4] + T[5]*S[5]; // check for corrupted tensor if (res>=0) { // estimate the bvalue from the base value and the norm of the gradient // - because of this estimation the vector have to be normalized beforehand // otherwise the modelled signal is wrong ( i.e. not scaled properly ) const double bval = m_BValue * twonorm * twonorm; out[i] = static_cast( maskvalue * 1.0 * b0 * exp ( -1.0 * bval * res ) ); } } } for(unsigned int idx = 0; idx < b0_indices.size(); idx++ ) { out[b0_indices.at(idx)] = b0; } itOut.Set(out); if( threadId==0 && step>0) { if( (progress%step)==0 ) { this->UpdateProgress ( double(progress)/double(numPixels) ); } } ++progress; ++itB0; ++itIn; ++itOut; } if( threadId==0 ) { this->UpdateProgress (1.0); } } +template +void +TensorImageToDiffusionImageFilter +::UpdateOutputInformation() +{ + // Calls to superclass updateoutputinformation + Superclass::UpdateOutputInformation(); + this->GetOutput()->SetVectorLength( m_GradientList->size() ); +} } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp index b4fb13c110..93dc5802af 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReader.cpp @@ -1,306 +1,311 @@ /*=================================================================== 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 "mitkDiffusionDICOMFileReader.h" #include "mitkDiffusionDICOMFileReaderHelper.h" #include "mitkDiffusionHeaderSiemensDICOMFileReader.h" #include "mitkDiffusionHeaderSiemensMosaicDICOMFileReader.h" #include "mitkDiffusionHeaderGEDICOMFileReader.h" #include "mitkDiffusionHeaderPhilipsDICOMFileReader.h" +#include +#include #include "mitkStringProperty.h" +#include static void PerformHeaderAnalysis( mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType headers ) { unsigned int images = headers.size(); unsigned int unweighted_images = 0; unsigned int weighted_images = 0; mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType::const_iterator c_iter = headers.begin(); while( c_iter != headers.end() ) { const mitk::DiffusionImageDICOMHeaderInformation h = *c_iter; if( h.baseline ) unweighted_images++; if( h.b_value > 0 ) weighted_images++; ++c_iter; } MITK_INFO << " :: Analyzed volumes " << images << "\n" << " :: \t"<< unweighted_images << " b = 0" << "\n" << " :: \t"<< weighted_images << " b > 0"; } mitk::DiffusionDICOMFileReader::DiffusionDICOMFileReader() { } mitk::DiffusionDICOMFileReader::~DiffusionDICOMFileReader() { } bool mitk::DiffusionDICOMFileReader ::LoadImages() { unsigned int numberOfOutputs = this->GetNumberOfOutputs(); bool success = true; for(unsigned int o = 0; o < numberOfOutputs; ++o) { success &= this->LoadSingleOutputImage( this->m_OutputHeaderContainer.at(o), this->InternalGetOutput(o), this->m_IsMosaicData.at(o) ); } return success; } bool mitk::DiffusionDICOMFileReader ::LoadSingleOutputImage( DiffusionHeaderDICOMFileReader::DICOMHeaderListType retrievedHeader, DICOMImageBlockDescriptor& block, bool is_mosaic) { // prepare data reading DiffusionDICOMFileReaderHelper helper; DiffusionDICOMFileReaderHelper::VolumeFileNamesContainer filenames; const DICOMImageFrameList& frames = block.GetImageFrameList(); int numberOfDWImages = block.GetIntProperty("timesteps", 1); int numberOfFramesPerDWImage = frames.size() / numberOfDWImages; assert( int( double((double) frames.size() / (double) numberOfDWImages)) == numberOfFramesPerDWImage ); for( int idx = 0; idx < numberOfDWImages; idx++ ) { std::vector< std::string > FileNamesPerVolume; DICOMImageFrameList::const_iterator timeStepStart = frames.begin() + idx * numberOfFramesPerDWImage; DICOMImageFrameList::const_iterator timeStepEnd = frames.begin() + (idx+1) * numberOfFramesPerDWImage; for (DICOMImageFrameList::const_iterator frameIter = timeStepStart; frameIter != timeStepEnd; ++frameIter) { FileNamesPerVolume.push_back( (*frameIter)->Filename ); } filenames.push_back( FileNamesPerVolume ); } // TODO : only prototyping to test loading of diffusion images // we need some solution for the different types - typedef mitk::DiffusionImage DiffusionImageType; - DiffusionImageType::Pointer output_image = DiffusionImageType::New(); + mitk::Image::Pointer output_image = mitk::Image::New(); - DiffusionImageType::GradientDirectionContainerType::Pointer directions = - DiffusionImageType::GradientDirectionContainerType::New(); + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer directions = + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::New(); double max_bvalue = 0; for( int idx = 0; idx < numberOfDWImages; idx++ ) { DiffusionImageDICOMHeaderInformation header = retrievedHeader.at(idx); if( max_bvalue < header.b_value ) max_bvalue = header.b_value; } // normalize the retrieved gradient directions according to the set b-value (maximal one) for( int idx = 0; idx < numberOfDWImages; idx++ ) { DiffusionImageDICOMHeaderInformation header = retrievedHeader.at(idx); - DiffusionImageType::GradientDirectionType grad = header.g_vector; + mitk::DiffusionPropertyHelper::GradientDirectionType grad = header.g_vector; grad.normalize(); grad *= sqrt( header.b_value / max_bvalue ); directions->push_back( grad ); } // initialize the output image - output_image->SetReferenceBValue( max_bvalue ); - output_image->SetDirections( directions ); + output_image->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( directions ) ); + output_image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( max_bvalue ) ); if( is_mosaic && this->m_ResolveMosaic ) { mitk::DiffusionHeaderSiemensMosaicDICOMFileReader::Pointer mosaic_reader = mitk::DiffusionHeaderSiemensMosaicDICOMFileReader::New(); // retrieve the remaining meta-information needed for mosaic reconstruction // it suffices to get it exemplatory from the first file in the file list mosaic_reader->RetrieveMosaicInformation( filenames.at(0).at(0) ); mitk::MosaicDescriptor mdesc = mosaic_reader->GetMosaicDescriptor(); - output_image->SetVectorImage( helper.LoadMosaicToVector( filenames, mdesc ) ); + mitk::CastToMitkImage( helper.LoadMosaicToVector( filenames, mdesc ), output_image ); } else { - output_image->SetVectorImage( helper.LoadToVector( filenames ) ); + mitk::CastToMitkImage( helper.LoadToVector( filenames ), output_image ); } - output_image->InitializeFromVectorImage(); + + mitk::DiffusionPropertyHelper propertyHelper( output_image ); + propertyHelper.InitializeImage(); + output_image->SetProperty("diffusion.dicom.importname", mitk::StringProperty::New( helper.GetOutputName(filenames) ) ); block.SetMitkImage( (mitk::Image::Pointer) output_image ); return block.GetMitkImage().IsNotNull(); } void mitk::DiffusionDICOMFileReader ::AnalyzeInputFiles() { this->SetGroup3DandT(true); Superclass::AnalyzeInputFiles(); // collect output from superclass size_t number_of_outputs = this->GetNumberOfOutputs(); if(number_of_outputs == 0) { MITK_ERROR << "Failed to parse input, retrieved 0 outputs from SeriesGDCMReader "; } MITK_INFO("diffusion.dicomreader") << "Retrieved " << number_of_outputs << "outputs."; for( unsigned int outputidx = 0; outputidx < this->GetNumberOfOutputs(); outputidx++ ) { DICOMImageBlockDescriptor block_0 = this->GetOutput(outputidx); // collect vendor ID from the first output, first image StringList inputFilename; DICOMImageFrameInfo::Pointer frame_0 = block_0.GetImageFrameList().at(0); inputFilename.push_back( frame_0->Filename ); mitk::DiffusionHeaderDICOMFileReader::Pointer headerReader; bool isMosaic = false; gdcm::Scanner gdcmScanner; gdcm::Tag t_vendor(0x008, 0x0070); gdcm::Tag t_imagetype(0x008, 0x008); // add DICOM Tag for vendor gdcmScanner.AddTag( t_vendor ); // add DICOM Tag for image type gdcmScanner.AddTag( t_imagetype ); if( gdcmScanner.Scan( inputFilename ) ) { // retrieve both vendor and image type std::string vendor = gdcmScanner.GetValue( frame_0->Filename.c_str(), t_vendor ); std::string image_type = gdcmScanner.GetValue( frame_0->Filename.c_str(), t_imagetype ); MITK_INFO("diffusion.dicomreader") << "Output " << outputidx+1 << " Got vendor: " << vendor << " image type " << image_type; // parse vendor tag if( vendor.find("SIEMENS") != std::string::npos ) //&& image_type.find("DIFFUSION") != std::string::npos ) { if( image_type.find("MOSAIC") != std::string::npos ) { headerReader = mitk::DiffusionHeaderSiemensMosaicDICOMFileReader::New(); isMosaic = true; } else { headerReader = mitk::DiffusionHeaderSiemensDICOMFileReader::New(); } } else if( vendor.find("GE") != std::string::npos ) { headerReader = mitk::DiffusionHeaderGEDICOMFileReader::New(); } else if( vendor.find("Philips") != std::string::npos ) { headerReader = mitk::DiffusionHeaderPhilipsDICOMFileReader::New(); } else { // unknown vendor } if( headerReader.IsNull() ) { MITK_ERROR << "No header reader for given vendor. "; continue; } } else { continue; } bool canread = false; // iterate over the threeD+t block int numberOfTimesteps = block_0.GetIntProperty("timesteps", 1); int framesPerTimestep = block_0.GetImageFrameList().size() / numberOfTimesteps; for( int idx = 0; idx < numberOfTimesteps; idx++ ) { int access_idx = idx * framesPerTimestep; DICOMImageFrameInfo::Pointer frame = this->GetOutput( outputidx ).GetImageFrameList().at( access_idx ); canread = headerReader->ReadDiffusionHeader( frame->Filename ); } if( canread ) { // collect the information mitk::DiffusionHeaderDICOMFileReader::DICOMHeaderListType retrievedHeader = headerReader->GetHeaderInformation(); m_IsMosaicData.push_back(isMosaic); m_OutputHeaderContainer.push_back( retrievedHeader ); m_OutputReaderContainer.push_back( headerReader ); } } this->SetNumberOfOutputs( this->m_OutputHeaderContainer.size() ); for( unsigned int outputidx = 0; outputidx < this->GetNumberOfOutputs(); outputidx++ ) { // TODO : Analyze outputs + header information, i.e. for the loading confidence MITK_INFO("diffusion.dicomreader") << "---- DICOM Analysis Report ---- :: Output " << outputidx+1 << " of " << this->GetNumberOfOutputs(); try{ PerformHeaderAnalysis( this->m_OutputHeaderContainer.at( outputidx) ); } catch( const std::exception& se) { MITK_ERROR << "STD Exception " << se.what(); } MITK_INFO("diffusion.dicomreader") << "==========================================="; } } bool mitk::DiffusionDICOMFileReader ::CanHandleFile(const std::string & /* filename */) { //FIXME : return true; } diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h index 6cd30e36d7..e4da4b1849 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionDICOMFileReaderHelper.h @@ -1,283 +1,281 @@ #ifndef MITKDIFFUSIONDICOMFILEREADERHELPER_H #define MITKDIFFUSIONDICOMFILEREADERHELPER_H -#include "mitkDiffusionImage.h" - #include "itkImageSeriesReader.h" #include "itkVectorImage.h" #include "itkImageRegionIteratorWithIndex.h" namespace mitk { /** * @brief The MosaicDescriptor struct is a help struct holding the necessary information for * loading a mosaic DICOM file into an MITK file with correct geometry information */ struct MosaicDescriptor { unsigned int nimages; bool slicenormalup; itk::ImageBase<3>::SpacingType spacing; itk::ImageBase<3>::DirectionType direction; float origin[3]; }; class DiffusionDICOMFileReaderHelper { public: typedef std::vector< std::string > StringContainer; typedef std::vector< StringContainer > VolumeFileNamesContainer; std::string GetOutputName(const VolumeFileNamesContainer& filenames) { typedef itk::Image< short, 3> InputImageType; typedef itk::ImageSeriesReader< InputImageType > SeriesReaderType; SeriesReaderType::Pointer probe_reader = SeriesReaderType::New(); probe_reader->SetFileNames( filenames.at(0) ); probe_reader->GenerateOutputInformation(); probe_reader->Update(); std::string seriesDescTag, seriesNumberTag, patientName; SeriesReaderType::DictionaryArrayRawPointer inputDict = probe_reader->GetMetaDataDictionaryArray(); if( ! itk::ExposeMetaData< std::string > ( *(*inputDict)[0], "0008|103e", seriesDescTag ) ) seriesDescTag = "UNSeries"; if( ! itk::ExposeMetaData< std::string > ( *(*inputDict)[0], "0020|0011", seriesNumberTag ) ) seriesNumberTag = "00000"; if( ! itk::ExposeMetaData< std::string > ( *(*inputDict)[0], "0010|0010", patientName ) ) patientName = "UnknownName"; std::stringstream ss; ss << seriesDescTag << "_" << seriesNumberTag << "_" << patientName; return ss.str(); } template< typename PixelType, unsigned int VecImageDimension> typename itk::VectorImage< PixelType, VecImageDimension >::Pointer LoadToVector( const VolumeFileNamesContainer& filenames //const itk::ImageBase<3>::RegionType requestedRegion ) { typedef itk::Image< PixelType, 3> InputImageType; typedef itk::ImageSeriesReader< InputImageType > SeriesReaderType; typename SeriesReaderType::Pointer probe_reader = SeriesReaderType::New(); probe_reader->SetFileNames( filenames.at(0) ); probe_reader->GenerateOutputInformation(); const itk::ImageBase<3>::RegionType requestedRegion = probe_reader->GetOutput()->GetLargestPossibleRegion(); MITK_INFO << " --- Probe reader : \n" << " Retrieved LPR " << requestedRegion; typedef itk::VectorImage< PixelType, 3 > VectorImageType; typename VectorImageType::Pointer output_image = VectorImageType::New(); output_image->SetNumberOfComponentsPerPixel( filenames.size() ); output_image->SetSpacing( probe_reader->GetOutput()->GetSpacing() ); output_image->SetOrigin( probe_reader->GetOutput()->GetOrigin() ); output_image->SetDirection( probe_reader->GetOutput()->GetDirection() ); output_image->SetLargestPossibleRegion( probe_reader->GetOutput()->GetLargestPossibleRegion() ); output_image->SetBufferedRegion( requestedRegion ); output_image->Allocate(); itk::ImageRegionIterator< VectorImageType > vecIter( output_image, requestedRegion ); VolumeFileNamesContainer::const_iterator volumesFileNamesIter = filenames.begin(); // iterate over the given volumes unsigned int component = 0; while( volumesFileNamesIter != filenames.end() ) { MITK_INFO << " ======== Loading volume " << component+1 << " of " << filenames.size(); typename SeriesReaderType::Pointer volume_reader = SeriesReaderType::New(); volume_reader->SetFileNames( *volumesFileNamesIter ); try { volume_reader->UpdateLargestPossibleRegion(); } catch( const itk::ExceptionObject &e) { mitkThrow() << " ITK Series reader failed : "<< e.what(); } itk::ImageRegionConstIterator< InputImageType > iRCIter ( volume_reader->GetOutput(), volume_reader->GetOutput()->GetLargestPossibleRegion() ); // transfer to vector image iRCIter.GoToBegin(); vecIter.GoToBegin(); while( !iRCIter.IsAtEnd() ) { typename VectorImageType::PixelType vector_pixel = vecIter.Get(); vector_pixel.SetElement( component, iRCIter.Get() ); vecIter.Set( vector_pixel ); ++vecIter; ++iRCIter; } ++volumesFileNamesIter; component++; } return output_image; } /** * Create the vector image for the resulting diffusion image from a mosaic DICOM image set, * The method needs to be provided with the MosaicDescriptor struct to be able to compute the * correct index and to set the geometry information of the image itself. */ template< typename PixelType, unsigned int VecImageDimension> typename itk::VectorImage< PixelType, VecImageDimension >::Pointer LoadMosaicToVector( const VolumeFileNamesContainer& filenames, const MosaicDescriptor& mosaicInfo ) { typedef itk::Image< PixelType, 3> MosaicImageType; typedef itk::ImageFileReader< MosaicImageType > SingleImageReaderType; // generate output typedef itk::VectorImage< PixelType, 3 > VectorImageType; VolumeFileNamesContainer::const_iterator volumesFileNamesIter = filenames.begin(); // probe the first file to retrieve the size of the 2d image // we need this information to compute the index relation between mosaic and resulting 3d position // but we need it only once typename SingleImageReaderType::Pointer mosaic_probe = SingleImageReaderType::New(); mosaic_probe->SetFileName( (*volumesFileNamesIter).at(0) ); try { mosaic_probe->UpdateLargestPossibleRegion(); } catch( const itk::ExceptionObject &e) { mitkThrow() << " ITK Image file reader failed : "<< e.what(); } typename MosaicImageType::RegionType mosaic_lpr = mosaic_probe->GetOutput()->GetLargestPossibleRegion(); MITK_INFO << " == MOSAIC: " << mosaic_lpr; itk::ImageBase<3>::SizeValueType images_per_row = ceil( sqrt( (float) mosaicInfo.nimages ) ); itk::ImageBase<3>::RegionType requestedRegion; requestedRegion.SetSize( 0, mosaic_lpr.GetSize()[0]/images_per_row); requestedRegion.SetSize( 1, mosaic_lpr.GetSize()[1]/images_per_row); requestedRegion.SetSize( 2, mosaicInfo.nimages); typename VectorImageType::Pointer output_image = VectorImageType::New(); output_image->SetNumberOfComponentsPerPixel( filenames.size() ); typename VectorImageType::DirectionType dmatrix; dmatrix.SetIdentity(); std::vector dirx = mosaic_probe->GetImageIO()->GetDirection(0); std::vector diry = mosaic_probe->GetImageIO()->GetDirection(1); std::vector dirz = mosaic_probe->GetImageIO()->GetDirection(2); dmatrix.GetVnlMatrix().set_column( 0, &dirx[0] ); dmatrix.GetVnlMatrix().set_column( 1, &diry[0] ); dmatrix.GetVnlMatrix().set_column( 2, &dirz[0] ); /* FIXME!!! The struct currently does not provide the geometry information the loading works as required*/ output_image->SetSpacing( mosaicInfo.spacing ); output_image->SetOrigin( mosaic_probe->GetOutput()->GetOrigin() ); output_image->SetDirection( dmatrix ); output_image->SetLargestPossibleRegion( requestedRegion ); output_image->SetBufferedRegion( requestedRegion ); output_image->Allocate(); itk::ImageRegionIteratorWithIndex< VectorImageType > vecIter( output_image, requestedRegion ); // hold the image sizes in an extra variable ( used very often ) typename MosaicImageType::SizeValueType dx = requestedRegion.GetSize()[0]; typename MosaicImageType::SizeValueType dy = requestedRegion.GetSize()[1]; // iterate over the given volumes unsigned int component = 0; while( volumesFileNamesIter != filenames.end() ) { MITK_INFO << " ======== Loading volume " << component+1 << " of " << filenames.size(); typename SingleImageReaderType::Pointer mosaic_reader = SingleImageReaderType::New(); mosaic_reader->SetFileName( (*volumesFileNamesIter).at(0) ); try { mosaic_reader->UpdateLargestPossibleRegion(); } catch( const itk::ExceptionObject &e) { mitkThrow() << " ITK Image file reader failed : "<< e.what(); } typename MosaicImageType::Pointer current_mosaic = mosaic_reader->GetOutput(); vecIter.GoToBegin(); while( !vecIter.IsAtEnd() ) { typename VectorImageType::PixelType vector_pixel = vecIter.Get(); typename VectorImageType::IndexType threeD_index = vecIter.GetIndex(); typename MosaicImageType::IndexType mosaic_index; mosaic_index[2] = 0; // first find the corresponding tile in the mosaic // this is defined by the z-position of the vector (3D) image iterator // in x : z_index % #images_in_grid // in y : z_index / #images_in_grid // // the remaining is just computing the correct position in the mosaic, done by // --------- index of (0,0,z) ----- + --- current 2d position --- mosaic_index[0] = (threeD_index[2] % images_per_row) * dx + threeD_index[0]; mosaic_index[1] = (threeD_index[2] / images_per_row) * dy + threeD_index[1]; typename MosaicImageType::PixelType mosaic_pixel = current_mosaic->GetPixel( mosaic_index ); vector_pixel.SetElement( component, mosaic_pixel ); vecIter.Set( vector_pixel ); ++vecIter; } ++volumesFileNamesIter; component++; } return output_image; } }; } #endif // MITKDIFFUSIONDICOMFILEREADERHELPER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h index ed62cf3f1e..7e1f348f36 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderDICOMFileReader.h @@ -1,117 +1,120 @@ /*=================================================================== 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 MITKDIFFUSIONHEADERFILEREADER_H #define MITKDIFFUSIONHEADERFILEREADER_H #include +#include #include -#include "mitkDiffusionImage.h" +#include #include "gdcmScanner.h" #include "gdcmReader.h" +#include + namespace mitk { /** * @brief The DiffusionImageHeaderInformation struct */ struct DiffusionImageDICOMHeaderInformation { DiffusionImageDICOMHeaderInformation() : b_value(0), baseline(false), isotropic(false) { g_vector.fill(0); } void Print() { MITK_INFO << " DiffusionImageHeaderInformation : \n" << " : b value : " << b_value << "\n" << " : gradient : " << g_vector << "\n" << " : isotropic : " << isotropic << "\n --- \n"; } unsigned int b_value; vnl_vector_fixed< double, 3> g_vector; bool baseline; bool isotropic; }; struct DiffusionImageMosaicDICOMHeaderInformation : public DiffusionImageDICOMHeaderInformation { unsigned long n_images; bool slicenormalup; }; /** * @class DiffusionHeaderDICOMFileReader * * @brief Abstract class for all vendor specific diffusion file header reader * * To provide a diffusion header reader for a new vendor, reimplement the \sa ReadDiffusionHeader method. */ class MitkDiffusionCore_EXPORT DiffusionHeaderDICOMFileReader : public itk::LightObject { public: typedef std::vector< DiffusionImageDICOMHeaderInformation > DICOMHeaderListType; mitkClassMacro( DiffusionHeaderDICOMFileReader, itk::LightObject ) itkSimpleNewMacro( Self ) /** * @brief IsDiffusionHeader Parse the given dicom file and collect the special diffusion image information * @return */ virtual bool ReadDiffusionHeader( std::string ){ return false; } DICOMHeaderListType GetHeaderInformation(); protected: DiffusionHeaderDICOMFileReader(); virtual ~DiffusionHeaderDICOMFileReader(); DICOMHeaderListType m_HeaderInformationList; }; /** * @brief Retrieve the value of a gdcm tag to the given string * * @param tag the gdcm::Tag to be search for * @param dataset a gdcm::DataSet to look into * @param target a string to store the value of the given tag if found * @param verbose make some output * * @return true if a string was found, false otherwise */ bool RevealBinaryTag(const gdcm::Tag tag, const gdcm::DataSet& dataset, std::string& target); bool RevealBinaryTagC(const gdcm::Tag tag, const gdcm::DataSet& dataset, char* target_array ); } // end namespace mitk #endif // MITKDIFFUSIONHEADERFILEREADER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp index 04f552619e..e90e9b77a4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp @@ -1,101 +1,102 @@ #include "mitkDiffusionHeaderGEDICOMFileReader.h" +#include mitk::DiffusionHeaderGEDICOMFileReader ::DiffusionHeaderGEDICOMFileReader() { } mitk::DiffusionHeaderGEDICOMFileReader ::~DiffusionHeaderGEDICOMFileReader() { } bool mitk::DiffusionHeaderGEDICOMFileReader ::ReadDiffusionHeader(std::string filename) { gdcm::Reader gdcmReader; gdcmReader.SetFileName( filename.c_str() ); gdcmReader.Read(); gdcm::Tag ge_bvalue_tag( 0x0043, 0x1039 ); gdcm::Tag ge_gradient_x( 0x0019, 0x10bb ); gdcm::Tag ge_gradient_y( 0x0019, 0x10bc ); gdcm::Tag ge_gradient_z( 0x0019, 0x10bd ); bool success = true; DiffusionImageDICOMHeaderInformation header_info; std::string ge_tagvalue_string; char* pEnd; // start with b-value success = RevealBinaryTag( ge_bvalue_tag, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string ); // b value stored in the first bytes // typical example: "1000\8\0\0" for bvalue=1000 // "40\8\0\0" for bvalue=40 // so we need to cut off the last 6 elements const char* bval_string = ge_tagvalue_string.substr(0,ge_tagvalue_string.length()-6).c_str(); header_info.b_value = static_cast(strtod( bval_string, &pEnd )); // now retrieve the gradient direction if(success && RevealBinaryTag( ge_gradient_x, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string ) ) { header_info.g_vector[0] = strtod( ge_tagvalue_string.c_str(), &pEnd ); } else { success = false; } if( success && RevealBinaryTag( ge_gradient_y, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string ) ) { header_info.g_vector[1] = strtod( ge_tagvalue_string.c_str(), &pEnd ); } else { success = false; } if( success && RevealBinaryTag( ge_gradient_z, gdcmReader.GetFile().GetDataSet(), ge_tagvalue_string ) ) { header_info.g_vector[2] = strtod( ge_tagvalue_string.c_str(), &pEnd ); } else { success = false; } if( success ) { // Fix for (0,0,0) direction in IVIM datasets if( header_info.b_value > 0 && header_info.g_vector.two_norm() < vnl_math::eps ) { header_info.g_vector.fill(1); header_info.g_vector.normalize(); header_info.isotropic = true; } // mark baseline if( header_info.b_value == 0 ) header_info.baseline = true; this->m_HeaderInformationList.push_back( header_info ); header_info.Print(); } return success; } diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp index b037bbb3da..a9c8710521 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp @@ -1,95 +1,97 @@ /*=================================================================== 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 "mitkDiffusionHeaderPhilipsDICOMFileReader.h" #include +#include + mitk::DiffusionHeaderPhilipsDICOMFileReader::DiffusionHeaderPhilipsDICOMFileReader() { } mitk::DiffusionHeaderPhilipsDICOMFileReader::~DiffusionHeaderPhilipsDICOMFileReader() { } bool mitk::DiffusionHeaderPhilipsDICOMFileReader::ReadDiffusionHeader(std::string filename) { gdcm::Reader gdcmReader; gdcmReader.SetFileName( filename.c_str() ); gdcmReader.Read(); gdcm::Tag philips_bvalue_tag( 0x2001, 0x1003 ); //gdcm::Tag philips_gradient_direction( 0x2001, 0x1004 ); DiffusionImageDICOMHeaderInformation header_info; //std::string tagvalue_string; //char* pEnd; // reveal b-value float bvalue = 0; if( RevealBinaryTagC( philips_bvalue_tag, gdcmReader.GetFile().GetDataSet(), (char*) &bvalue) ) { header_info.b_value = std::ceil( bvalue ); if( header_info.b_value == 0) header_info.baseline = true; } else { MITK_WARN("diffusion.dicomreader.philips") << "No b-value found. Most probably no diffusion-weighted image."; return false; } gdcm::Tag philips_new_bvalue_tag( 0x0018,0x9087 ); double dbvalue = 0; if( RevealBinaryTagC( philips_new_bvalue_tag, gdcmReader.GetFile().GetDataSet(), (char*) &dbvalue) ) { MITK_INFO("philips.dicom.diffusion.bvalue") << dbvalue; } if( header_info.baseline ) { // no direction in unweighted images header_info.g_vector.fill(0); } else { MITK_INFO("philips.dicom.diffusion.gradientdir") << "Parsing gradient direction."; gdcm::Tag philips_gradient_direction_new( 0x0018, 0x9089 ); double gr_dir_arr[3] = {1,0,-1}; if( RevealBinaryTagC( philips_gradient_direction_new, gdcmReader.GetFile().GetDataSet(), (char*) &gr_dir_arr ) ) { MITK_INFO("philips.dicom.diffusion.gradient") << "(" << gr_dir_arr[0] <<"," << gr_dir_arr[1] <<"," <m_HeaderInformationList.push_back( header_info ); return true; } diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp index 043f5e94f2..d5e17f2e14 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp @@ -1,42 +1,42 @@ #include "mitkDiffusionHeaderSiemensDICOMFileHelper.h" #include #include mitk::SiemensDiffusionHeaderType mitk::GetHeaderType( std::string header ) { // The CSA2 format begins with the string ‘SV10’, the CSA1 format does not. if( header.find("SV10") != std::string::npos ) { return mitk::SIEMENS_CSA2; } else { return mitk::SIEMENS_CSA1; } } bool mitk::ParseInputString( std::string input, std::vector& values, Siemens_Header_Format format_specs ) { - // TODO : Compute offset based on the format_specs, where does the 84 come from??? int offset = 84; int vm = *(input.c_str() + format_specs.NameLength ); for (int k = 0; k < vm; k++) { int itemLength = *(input.c_str() + offset + 4); int strideSize = static_cast (ceil(static_cast(itemLength)/4) * 4); std::string valueString = input.substr( offset+16, itemLength ); double value = atof( valueString.c_str() ); values.push_back( value ); offset += 16+strideSize; } - return true; + // If there are no values it is invalid + return (values.size() > 0 ); } diff --git a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp index 62416888ca..5cf26e5b8e 100644 --- a/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp @@ -1,151 +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. ===================================================================*/ #include "mitkDiffusionHeaderSiemensDICOMFileReader.h" #include "mitkDiffusionHeaderSiemensDICOMFileHelper.h" #include "gdcmScanner.h" #include "gdcmReader.h" +#include /** * @brief Extract b value from the siemens diffusion tag */ bool mitk::DiffusionHeaderSiemensDICOMFileReader ::ExtractSiemensDiffusionTagInformation( std::string tag_value, mitk::DiffusionImageDICOMHeaderInformation& values) { SiemensDiffusionHeaderType hformat = mitk::GetHeaderType( tag_value ); Siemens_Header_Format specs = this->m_SiemensFormatsCollection.at( hformat ); MITK_DEBUG << " Header format: " << hformat; MITK_DEBUG << " :: Retrieving b value. "; std::string::size_type tag_position = tag_value.find( "B_value", 0 ); if( tag_position == std::string::npos ) { MITK_ERROR << "No b value information found. "; return false; } std::string value_string = tag_value.substr( tag_position, tag_value.size() - tag_position + 1 ); std::vector value_array; if( ParseInputString(value_string, value_array, specs) ) { values.b_value = value_array.at(0); } else { MITK_INFO("diffusion.dicomreader.siemens") << "No b-value tag found. "; return false; } // search for GradientDirectionInformation if the bvalue is not null if( values.b_value > 0 ) { std::string::size_type tag_position = tag_value.find( "DiffusionGradientDirection", 0 ); // Possibly it is a IVIM dataset, i.e. the gradient direction is not relevant // and possibly either not set or set to zero if( tag_position == std::string::npos ) { MITK_WARN << "No gradient direction information, but non-zero b-value. Possibly an IVIM dataset. " << "\n" << "Setting gradient to (1,1,1)."; values.isotropic = true; values.g_vector.fill(1); return false; } value_array.clear(); std::string gradient_direction_str = tag_value.substr( tag_position, tag_value.size() - tag_position + 1 ); if( ParseInputString(gradient_direction_str, value_array, specs) ) { if( value_array.size() != 3 ) { MITK_ERROR << " Retrieved gradient information of length " << value_array.size(); return false; } for( unsigned int i=0; iExtractSiemensDiffusionTagInformation( siemens_diffusionheader_str, values ); m_HeaderInformationList.push_back( values ); } return true; } diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp index 2f909c5b73..d93af7cf2d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp @@ -1,117 +1,118 @@ /*=================================================================== 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 -template< typename TPixelType > -mitk::DiffusionImageCorrectionFilter::DiffusionImageCorrectionFilter() + +mitk::DiffusionImageCorrectionFilter::DiffusionImageCorrectionFilter() { } -template< typename TPixelType > -typename mitk::DiffusionImageCorrectionFilter::TransformMatrixType -mitk::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; } -template< typename TPixelType > -void mitk::DiffusionImageCorrectionFilter + +void mitk::DiffusionImageCorrectionFilter ::CorrectDirections( const TransformsVectorType& transformations) { if( m_SourceImage.IsNull() ) { mitkThrow() << " No diffusion image given! "; } - GradientDirectionContainerPointerType directions = m_SourceImage->GetDirections(); + 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->SetDirections( corrected_directions ); + m_SourceImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( corrected_directions ) ); } -template< typename TPixelType > -void mitk::DiffusionImageCorrectionFilter + +void mitk::DiffusionImageCorrectionFilter ::CorrectDirections( const TransformMatrixType& transformation) { if( m_SourceImage.IsNull() ) { mitkThrow() << " No diffusion image given! "; } TransformsVectorType transfVec; - for (unsigned int i=0; i< m_SourceImage->GetDirections()->Size();i++) + 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/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h index 7d6c3b3e59..405cea700d 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.h @@ -1,97 +1,96 @@ /*=================================================================== 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 "mitkDiffusionImageSource.h" +#include "mitkImageSource.h" +#include namespace mitk { /** * @class DiffusionImageCorrectionFilter */ -template< typename TPixelType > -class DiffusionImageCorrectionFilter - : public DiffusionImageSource< TPixelType > +class MitkDiffusionCore_EXPORT DiffusionImageCorrectionFilter + : public ImageSource { public: /** class macros */ mitkClassMacro( DiffusionImageCorrectionFilter, - DiffusionImageSource ) + 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 typename Superclass::OutputType DiffusionImageType; - typedef typename DiffusionImageType::Pointer DiffusionImageTypePointer; - typedef itk::VectorImage ImageType; + 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& ); 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; }; } -#include "mitkDiffusionImageCorrectionFilter.cpp" - #endif // MITKDIFFUSIONIMAGECORRECTIONFILTER_H diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake b/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake index 4a2a9d22b0..9e87b8b605 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/files.cmake @@ -1,16 +1,15 @@ set(MODULE_TESTS - mitkDiffusionImageEqualTest.cpp mitkNonLocalMeansDenoisingTest.cpp - mitkDiffusionImageEqualTest.cpp + mitkDiffusionPropertySerializerTest.cpp ) set(MODULE_CUSTOM_TESTS mitkPyramidImageRegistrationMethodTest.cpp mitkDWHeadMotionCorrectionTest.cpp mitkImageReconstructionTest.cpp mitkConvertDWITypeTest.cpp mitkExtractSingleShellTest.cpp mitkNonLocalMeansDenoisingTest.cpp mitkDiffusionDICOMFileReaderTest.cpp ) diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp index c157efef62..101143f9d9 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkB0ExtractionToSeparateImagesFilterTest.cpp @@ -1,109 +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 "itkVectorImage.h" #include "mitkBaseData.h" #include "mitkBaseDataIOFactory.h" #include "mitkImageWriter.h" #include "itkImageFileWriter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" -#include "mitkDiffusionImage.h" #include "mitkDiffusionCoreObjectFactory.h" /** Documentation * Test for factory registration */ int mitkB0ExtractionToSeparateImagesFilterTest(int argc , char* argv[]) { // always start with this! MITK_TEST_BEGIN("mitkB0ExtractionToSeparateImagesFilterTest"); MITK_TEST_CONDITION_REQUIRED(argc > 2, "Test image is specified. "); typedef short DiffusionPixelType; - typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; + typedef mitk::Image DiffusionImageType; typedef mitk::NrrdDiffusionImageReader< DiffusionPixelType> DiffusionNrrdReaderType; RegisterDiffusionImagingObjectFactory(); std::string inputFileName( argv[1] ); std::vector inputBaseDataVector = mitk::BaseDataIO::LoadBaseDataFromFile( inputFileName, "","",false); MITK_TEST_CONDITION_REQUIRED( inputBaseDataVector.size() > 0, "BaseDataIO returned non-empty vector."); mitk::BaseData::Pointer baseData = inputBaseDataVector.at(0); MITK_TEST_CONDITION_REQUIRED( baseData.IsNotNull(), "BaseData is not null") DiffusionImageType* vols = dynamic_cast< DiffusionImageType* >(baseData.GetPointer()); MITK_TEST_CONDITION_REQUIRED( vols != NULL, "Casting basedata to diffusion image successfull." ); // filter typedef itk::B0ImageExtractionToSeparateImageFilter< short, short> FilterType; typename FilterType::Pointer filter = FilterType::New(); MITK_TEST_CONDITION_REQUIRED(filter.IsNotNull(), "Filter instance created. "); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); filter->Update(); // output mitk::Image::Pointer mitkImage = mitk::Image::New(); MITK_TEST_CONDITION_REQUIRED( mitkImage.IsNotNull(), "mitkImage not null." ); mitkImage->InitializeByItk( filter->GetOutput() ); MITK_TEST_CONDITION_REQUIRED( mitkImage->GetDimension()==4, "Output image is a 4D image."); mitkImage->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); typedef itk::ImageFileWriter< FilterType::OutputImageType > itkImageWriterType; typename itkImageWriterType::Pointer itkWriter = itkImageWriterType::New(); itkWriter->SetFileName( argv[2] ); itkWriter->SetInput( filter->GetOutput() ); try { itkWriter->Update(); } catch(itk::ExceptionObject &e) { MITK_ERROR << "Catched exception from image writer. " << e.what(); } /* // write output mitk::ImageWriter::Pointer writer = mitk::ImageWriter::New(); MITK_TEST_CONDITION_REQUIRED( writer.IsNotNull(), "Writer instance created. "); writer->SetInput( mitkImage ); writer->SetExtension(".nrrd"); writer->SetFileName( "/localdata/hering/_Images/TestB0Extraction" ); writer->Update(); */ // always end with this! MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp index aeae0f8536..4e6449576a 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkConvertDWITypeTest.cpp @@ -1,51 +1,46 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "mitkIOUtil.h" #include "mitkDWIHeadMotionCorrectionFilter.h" -typedef short DiffusionPixelType; -typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; - /** * @brief Custom test to provide CMD-line access to the mitk::DWIHeadMotionCorrectionFilter * * @param argv : Input and Output image full path */ int mitkConvertDWITypeTest( int argc, char* argv[] ) { MITK_TEST_BEGIN("mitkConvertDWITypeTest"); MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output."); mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); - DiffusionImageType* dwimage = - static_cast( inputImage.GetPointer() ); try { - mitk::IOUtil::Save(dwimage, argv[2]); + mitk::IOUtil::Save(inputImage, argv[2]); } catch( const itk::ExceptionObject& e) { - MITK_ERROR << "Catched exception: " << e.what(); + MITK_ERROR << "Caught 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 7805c72b3c..e19fc61f9c 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDWHeadMotionCorrectionTest.cpp @@ -1,60 +1,60 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "mitkIOUtil.h" -#include "mitkDWIHeadMotionCorrectionFilter.h" +#include typedef short DiffusionPixelType; -typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; +typedef mitk::Image DiffusionImageType; /** * @brief Custom test to provide CMD-line access to the mitk::DWIHeadMotionCorrectionFilter * * @param argv : Input and Output image full path */ int mitkDWHeadMotionCorrectionTest( int argc, char* argv[] ) { MITK_TEST_BEGIN("mitkDWHeadMotionCorrectionTest"); MITK_TEST_CONDITION_REQUIRED( argc > 2, "Specify input and output."); // itk::MultiThreader::SetGlobalMaximumNumberOfThreads(1); mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); DiffusionImageType* dwimage = static_cast( inputImage.GetPointer() ); - mitk::DWIHeadMotionCorrectionFilter::Pointer corrfilter = - mitk::DWIHeadMotionCorrectionFilter::New(); + mitk::DWIHeadMotionCorrectionFilter::Pointer corrfilter = + mitk::DWIHeadMotionCorrectionFilter::New(); corrfilter->SetInput( dwimage ); corrfilter->SetAverageUnweighted(false); corrfilter->Update(); try { mitk::IOUtil::SaveBaseData(corrfilter->GetOutput(), argv[2]); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Catched exception: " << e.what(); mitkThrow() << "Failed with exception from subprocess!"; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp index 1b18949742..bfdd73c944 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionDICOMFileReaderTest.cpp @@ -1,115 +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. ===================================================================*/ #include "mitkDiffusionDICOMFileReader.h" #include "mitkDiffusionDICOMFileReaderTestHelper.h" #include "mitkDICOMTagBasedSorter.h" #include "mitkDICOMSortByTag.h" #include #include "mitkTestingMacros.h" using mitk::DICOMTag; int mitkDiffusionDICOMFileReaderTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkDiffusionDICOMFileReaderTest"); mitk::DiffusionDICOMFileReader::Pointer gdcmReader = mitk::DiffusionDICOMFileReader::New(); MITK_TEST_CONDITION_REQUIRED(gdcmReader.IsNotNull(), "DICOMITKSeriesGDCMReader can be instantiated."); std::string output_filename = "/tmp/dicom_out.dwi"; if( argc > 3) { mitk::DICOMFileReaderTestHelper::SetTestInputFilenames( argc-1,argv ); output_filename = std::string( argv[argc-1] ); } else { mitk::DICOMFileReaderTestHelper::SetTestInputFilenames( argc,argv ); } // check the Set/GetInput function mitk::DICOMFileReaderTestHelper::TestInputFilenames( gdcmReader ); MITK_INFO << "Test input filenanems"; // check that output is a good reproduction of input (no duplicates, no new elements) mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); MITK_INFO << "Test output"; // repeat test with some more realistic sorting gdcmReader = mitk::DiffusionDICOMFileReader::New(); // this also tests destruction mitk::DICOMTagBasedSorter::Pointer tagSorter = mitk::DICOMTagBasedSorter::New(); // Use tags as in Qmitk // all the things that split by tag in DicomSeriesReader tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0010) ); // Number of Rows tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0011) ); // Number of Columns tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0030) ); // Pixel Spacing tagSorter->AddDistinguishingTag( DICOMTag(0x0018, 0x1164) ); // Imager Pixel Spacing tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0037) ); // Image Orientation (Patient) // TODO add tolerance parameter (l. 1572 of original code) // TODO handle as real vectors! cluster with configurable errors! //tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x000e) ); // Series Instance UID //tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0010) ); tagSorter->AddDistinguishingTag( DICOMTag(0x0018, 0x0050) ); // Slice Thickness tagSorter->AddDistinguishingTag( DICOMTag(0x0028, 0x0008) ); // Number of Frames tagSorter->AddDistinguishingTag( DICOMTag(0x0020, 0x0052) ); // Frame of Reference UID // gdcmReader->AddSortingElement( tagSorter ); //mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); mitk::DICOMSortCriterion::ConstPointer sorting = mitk::DICOMSortByTag::New( DICOMTag(0x0020, 0x0013), // instance number mitk::DICOMSortByTag::New( DICOMTag(0x0020, 0x0012), // aqcuisition number mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0032), // aqcuisition time mitk::DICOMSortByTag::New( DICOMTag(0x0018, 0x1060), // trigger time mitk::DICOMSortByTag::New( DICOMTag(0x0008, 0x0018) // SOP instance UID (last resort, not really meaningful but decides clearly) ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer() ).GetPointer(); tagSorter->SetSortCriterion( sorting ); MITK_INFO << "Created sort"; gdcmReader->AddSortingElement( tagSorter ); mitk::DICOMFileReaderTestHelper::TestOutputsContainInputs( gdcmReader ); MITK_INFO << "Created sort"; //gdcmReader->PrintOutputs(std::cout, true); // really load images //mitk::DICOMFileReaderTestHelper::TestMitkImagesAreLoaded( gdcmReader ); gdcmReader->LoadImages(); mitk::Image::Pointer loaded_image = gdcmReader->GetOutput(0).GetMitkImage(); - mitk::DiffusionImage::Pointer d_img = static_cast*>( loaded_image.GetPointer() ); - try { - mitk::IOUtil::Save(d_img, output_filename.c_str()); + mitk::IOUtil::Save(loaded_image, output_filename.c_str()); } catch( const itk::ExceptionObject& e) { MITK_TEST_FAILED_MSG( << "Writer failed : " << e.what() ); } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp new file mode 100644 index 0000000000..436646692c --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkDiffusionPropertySerializerTest.cpp @@ -0,0 +1,184 @@ +/*=================================================================== + +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 "mitkIOUtil.h" +#include "mitkTestingMacros.h" +#include "mitkTestFixture.h" + +#include + +#include + +#include +#include + +#include + +class mitkDiffusionPropertySerializerTestSuite : public mitk::TestFixture +{ + + CPPUNIT_TEST_SUITE(mitkDiffusionPropertySerializerTestSuite); + MITK_TEST(Equal_SerializeandDeserialize_ReturnsTrue); + + //MITK_TEST(Equal_DifferentChannels_ReturnFalse); + + + CPPUNIT_TEST_SUITE_END(); + +private: + + /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ + mitk::PropertyList::Pointer propList; //represet image propertylist + mitk::BValueMapProperty::Pointer bvaluemap_prop; + mitk::GradientDirectionsProperty::Pointer gradientdirection_prop; + mitk::MeasurementFrameProperty::Pointer measurementframe_prop; + +public: + + /** +* @brief Setup Always call this method before each Test-case to ensure correct and new intialization of the used members for a new test case. (If the members are not used in a test, the method does not need to be called). +*/ + void setUp() + { + + propList = mitk::PropertyList::New(); + + mitk::BValueMapProperty::BValueMap map; + std::vector indices1; + indices1.push_back(1); + indices1.push_back(2); + indices1.push_back(3); + indices1.push_back(4); + + map[0] = indices1; + std::vector indices2; + indices2.push_back(4); + indices2.push_back(3); + indices2.push_back(2); + indices2.push_back(1); + + map[1000] = indices2; + bvaluemap_prop = mitk::BValueMapProperty::New(map).GetPointer(); + propList->SetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str(),bvaluemap_prop); + + + mitk::GradientDirectionsProperty::GradientDirectionsContainerType::Pointer gdc; + gdc = mitk::GradientDirectionsProperty::GradientDirectionsContainerType::New(); + + double a[3] = {3.0,4.0,1.4}; + vnl_vector_fixed vec1; + vec1.set(a); + gdc->push_back(vec1); + + double b[3] = {1.0,5.0,123.4}; + vnl_vector_fixed vec2; + vec2.set(b); + gdc->push_back(vec2); + + double c[3] = {13.0,84.02,13.4}; + vnl_vector_fixed vec3; + vec3.set(c); + gdc->push_back(vec3); + + + gradientdirection_prop = mitk::GradientDirectionsProperty::New(gdc.GetPointer()).GetPointer(); + propList->SetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), gradientdirection_prop); + + mitk::MeasurementFrameProperty::MeasurementFrameType mft; + + double row0[3] = {1,0,0}; + double row1[3] = {0,1,0}; + double row2[3] = {0,0,1}; + + mft.set_row(0,row0); + mft.set_row(1,row1); + mft.set_row(2,row2); + + measurementframe_prop = mitk::MeasurementFrameProperty::New(mft).GetPointer(); + propList->SetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), measurementframe_prop); + } + + void tearDown() + { + + } + + void Equal_SerializeandDeserialize_ReturnsTrue() + { + + assert(propList); + + /* try to serialize each property in the list, then deserialize again and check for equality */ + for (mitk::PropertyList::PropertyMap::const_iterator it = propList->GetMap()->begin(); it != propList->GetMap()->end(); ++it) + { + const mitk::BaseProperty* prop = it->second; + // construct name of serializer class + std::string serializername = std::string(prop->GetNameOfClass()) + "Serializer"; + std::list allSerializers = itk::ObjectFactoryBase::CreateAllInstance(serializername.c_str()); + MITK_TEST_CONDITION(allSerializers.size() > 0, std::string("Creating serializers for ") + serializername); + if (allSerializers.size() == 0) + { + MITK_TEST_OUTPUT( << "serialization not possible, skipping " << prop->GetNameOfClass()); + continue; + } + if (allSerializers.size() > 1) + { + MITK_TEST_OUTPUT (<< "Warning: " << allSerializers.size() << " serializers found for " << prop->GetNameOfClass() << "testing only the first one."); + } + mitk::BasePropertySerializer* serializer = dynamic_cast( allSerializers.begin()->GetPointer()); + MITK_TEST_CONDITION(serializer != NULL, serializername + std::string(" is valid")); + if (serializer != NULL) + { + serializer->SetProperty(prop); + TiXmlElement* valueelement = NULL; + try + { + valueelement = serializer->Serialize(); +// TiXmlPrinter p; +// valueelement->Accept(&p); +// MITK_INFO << p.CStr(); + } + catch (...) + { + } + MITK_TEST_CONDITION(valueelement != NULL, std::string("Serialize property with ") + serializername); + + if (valueelement == NULL) + { + MITK_TEST_OUTPUT( << "serialization failed, skipping deserialization"); + continue; + } + + mitk::BaseProperty::Pointer deserializedProp = serializer->Deserialize( valueelement ); + MITK_TEST_CONDITION(deserializedProp.IsNotNull(), "serializer created valid property"); + if (deserializedProp.IsNotNull()) + { + MITK_TEST_CONDITION(*(deserializedProp.GetPointer()) == *prop, "deserialized property equals initial property for type " << prop->GetNameOfClass()); + } + + } + else + { + MITK_TEST_OUTPUT( << "created serializer object is of class " << allSerializers.begin()->GetPointer()->GetNameOfClass()) + } + } // for all properties + + } + + +}; + +MITK_TEST_SUITE_REGISTRATION(mitkDiffusionPropertySerializer) diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp index f2b1ae2338..43686ea03e 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkExtractSingleShellTest.cpp @@ -1,96 +1,103 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkTestingMacros.h" #include "mitkIOUtil.h" #include #include "mitkDWIHeadMotionCorrectionFilter.h" +#include +#include +#include typedef short DiffusionPixelType; -typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; int mitkExtractSingleShellTest( int argc, char* argv[] ) { MITK_TEST_BEGIN("mitkExtractSingleShellTest"); MITK_TEST_CONDITION_REQUIRED( argc > 3, "Specify input and output and the shell to be extracted"); /* 1. Get input data */ - mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage( argv[1] ); - DiffusionImageType* dwimage = - static_cast( inputImage.GetPointer() ); + mitk::Image::Pointer dwimage = mitk::IOUtil::LoadImage( argv[1] ); - MITK_TEST_CONDITION_REQUIRED( dwimage != NULL, "Input is a dw-image"); + mitk::GradientDirectionsProperty::Pointer gradientsProperty = static_cast( dwimage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ); + + MITK_TEST_CONDITION_REQUIRED( gradientsProperty.IsNotNull(), "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; + typedef mitk::DiffusionPropertyHelper::BValueMapType BValueMap; // GetShellSelection from GUI BValueMap shellSelectionMap; - BValueMap originalShellMap = dwimage->GetBValueMap(); + BValueMap originalShellMap = static_cast(dwimage->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap(); std::vector newNumGradientDirections; shellSelectionMap[extract_value] = originalShellMap[extract_value]; newNumGradientDirections.push_back( originalShellMap[extract_value].size() ) ; - DiffusionImageType::GradientDirectionContainerType::Pointer gradientContainer = dwimage->GetDirections(); + itk::VectorImage< short, 3 >::Pointer itkVectorImagePointer = itk::VectorImage< short, 3 >::New(); + mitk::CastToItkImage(dwimage, itkVectorImagePointer); + itk::VectorImage< short, 3 > *vectorImage = itkVectorImagePointer.GetPointer(); + + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradientContainer = static_cast( dwimage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); FilterType::Pointer filter = FilterType::New(); - filter->SetInput(dwimage->GetVectorImage()); + filter->SetInput(vectorImage); 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->SetReferenceBValue(dwimage->GetReferenceBValue()); - image->SetDirections(filter->GetGradientDirections()); - image->SetMeasurementFrame(dwimage->GetMeasurementFrame()); - image->InitializeFromVectorImage(); + mitk::Image::Pointer outImage = mitk::GrabItkImageMemory( filter->GetOutput() ); + outImage->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( filter->GetGradientDirections() ) ); + outImage->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(dwimage->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); + outImage->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( static_cast(dwimage->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() )->GetMeasurementFrame() ) ); + mitk::DiffusionPropertyHelper propertyHelper( outImage ); + propertyHelper.InitializeImage(); + /* * 3. Write output data **/ try { - mitk::IOUtil::Save(image, argv[2]); + mitk::IOUtil::Save(outImage, argv[2]); } catch( const itk::ExceptionObject& e) { MITK_ERROR << "Catched exception: " << e.what(); mitkThrow() << "Failed with exception from subprocess!"; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp index e1f16dcca7..33bfd4cd53 100755 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkImageReconstructionTest.cpp @@ -1,152 +1,158 @@ /*=================================================================== 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 #include #include -#include #include #include #include #include #include #include #include #include +#include +#include using namespace std; int mitkImageReconstructionTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkImageReconstructionTest"); MITK_TEST_CONDITION_REQUIRED(argc>1,"check for input data") try { - mitk::DiffusionImage::Pointer dwi = dynamic_cast*>(mitk::IOUtil::LoadDataNode(argv[1])->GetData()); + mitk::Image::Pointer dwi = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[1])->GetData()); + itk::VectorImage::Pointer itkVectorImagePointer = itk::VectorImage::New(); + mitk::CastToItkImage(dwi, itkVectorImagePointer); + + float b_value = mitk::DiffusionPropertyHelper::GetReferenceBValue( dwi ); + mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer gradients = mitk::DiffusionPropertyHelper::GetGradientContainer(dwi); { MITK_INFO << "Tensor reconstruction " << argv[2]; mitk::TensorImage::Pointer tensorImage = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[2])->GetData()); typedef itk::DiffusionTensor3DReconstructionImageFilter< short, short, float > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer filter = TensorReconstructionImageFilterType::New(); - filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); - filter->SetBValue(dwi->GetReferenceBValue()); + filter->SetGradientImage( gradients, itkVectorImagePointer ); + filter->SetBValue( b_value ); filter->Update(); mitk::TensorImage::Pointer testImage = mitk::TensorImage::New(); testImage->InitializeByItk( filter->GetOutput() ); testImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *tensorImage, 0.0001, true), "tensor reconstruction test."); } { MITK_INFO << "Numerical Q-ball reconstruction " << argv[3]; mitk::QBallImage::Pointer qballImage = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[3])->GetData()); typedef itk::DiffusionQballReconstructionImageFilter QballReconstructionImageFilterType; QballReconstructionImageFilterType::Pointer filter = QballReconstructionImageFilterType::New(); - filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); - filter->SetBValue(dwi->GetReferenceBValue()); + filter->SetGradientImage( gradients, itkVectorImagePointer ); + filter->SetBValue( b_value ); filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); filter->Update(); mitk::QBallImage::Pointer testImage = mitk::QBallImage::New(); testImage->InitializeByItk( filter->GetOutput() ); testImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "Numerical Q-ball reconstruction test."); } { MITK_INFO << "Standard Q-ball reconstruction " << argv[4]; mitk::QBallImage::Pointer qballImage = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[4])->GetData()); typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); - filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); - filter->SetBValue(dwi->GetReferenceBValue()); + filter->SetGradientImage( gradients, itkVectorImagePointer ); + filter->SetBValue( b_value ); filter->SetLambda(0.006); filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); filter->Update(); mitk::QBallImage::Pointer testImage = mitk::QBallImage::New(); testImage->InitializeByItk( filter->GetOutput() ); testImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "Standard Q-ball reconstruction test."); } { MITK_INFO << "CSA Q-ball reconstruction " << argv[5]; mitk::QBallImage::Pointer qballImage = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[5])->GetData()); typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); - filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); - filter->SetBValue(dwi->GetReferenceBValue()); + filter->SetGradientImage( gradients, itkVectorImagePointer ); + filter->SetBValue( b_value ); filter->SetLambda(0.006); filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); filter->Update(); mitk::QBallImage::Pointer testImage = mitk::QBallImage::New(); testImage->InitializeByItk( filter->GetOutput() ); testImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "CSA Q-ball reconstruction test."); } { MITK_INFO << "ADC profile reconstruction " << argv[6]; mitk::QBallImage::Pointer qballImage = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[6])->GetData()); typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); - filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); - filter->SetBValue(dwi->GetReferenceBValue()); + filter->SetGradientImage( gradients, itkVectorImagePointer ); + filter->SetBValue( b_value ); filter->SetLambda(0.006); filter->SetNormalizationMethod(FilterType::QBAR_ADC_ONLY); filter->Update(); mitk::QBallImage::Pointer testImage = mitk::QBallImage::New(); testImage->InitializeByItk( filter->GetOutput() ); testImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "ADC profile reconstruction test."); } { MITK_INFO << "Raw signal modeling " << argv[7]; mitk::QBallImage::Pointer qballImage = dynamic_cast(mitk::IOUtil::LoadDataNode(argv[7])->GetData()); typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); - filter->SetGradientImage( dwi->GetDirections(), dwi->GetVectorImage() ); - filter->SetBValue(dwi->GetReferenceBValue()); + filter->SetGradientImage( gradients, itkVectorImagePointer ); + filter->SetBValue( b_value ); filter->SetLambda(0.006); filter->SetNormalizationMethod(FilterType::QBAR_RAW_SIGNAL); filter->Update(); mitk::QBallImage::Pointer testImage = mitk::QBallImage::New(); testImage->InitializeByItk( filter->GetOutput() ); testImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); MITK_TEST_CONDITION_REQUIRED(mitk::Equal(*testImage, *qballImage, 0.0001, true), "Raw signal modeling test."); } } catch (itk::ExceptionObject e) { MITK_INFO << e; return EXIT_FAILURE; } catch (std::exception e) { MITK_INFO << e.what(); return EXIT_FAILURE; } catch (...) { MITK_INFO << "ERROR!?!"; return EXIT_FAILURE; } MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp index 680b5c6295..62a7941cb4 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp @@ -1,176 +1,169 @@ /*=================================================================== 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 "mitkDiffusionImage.h" #include "mitkIOUtil.h" #include "mitkTestingMacros.h" #include "mitkTestFixture.h" #include "itkNonLocalMeansDenoisingFilter.h" +#include "mitkGradientDirectionsProperty.h" +#include "mitkITKImageImport.h" +#include class mitkNonLocalMeansDenoisingTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkNonLocalMeansDenoisingTestSuite); MITK_TEST(Denoise_NLMg_shouldReturnTrue); MITK_TEST(Denoise_NLMr_shouldReturnTrue); MITK_TEST(Denoise_NLMv_shouldReturnTrue); MITK_TEST(Denoise_NLMvr_shouldReturnTrue); CPPUNIT_TEST_SUITE_END(); private: + typedef itk::VectorImage VectorImagetType; + /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ - mitk::DiffusionImage::Pointer m_Image; - mitk::DiffusionImage::Pointer m_ReferenceImage; - mitk::DiffusionImage::Pointer m_DenoisedImage; + mitk::Image::Pointer m_Image; + mitk::Image::Pointer m_ReferenceImage; + mitk::Image::Pointer m_DenoisedImage; itk::Image::Pointer m_ImageMask; itk::NonLocalMeansDenoisingFilter::Pointer m_DenoisingFilter; public: /** * @brief Setup Always call this method before each Test-case to ensure correct and new intialization of the used members for a new test case. (If the members are not used in a test, the method does not need to be called). */ void setUp() { //generate test images std::string imagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi.dwi"); -// std::string maskPath = GetTestDataFilePath("DiffusionImaging/Denoising/denoising_mask.nrrd"); - m_Image = static_cast*>( mitk::IOUtil::LoadImage(imagePath).GetPointer()); -// mitk::Image::Pointer imageMask = static_cast( mitk::IOUtil::LoadImage(maskPath).GetPointer()); -// mitk::CastToItkImage(imageMask, m_ImageMask); + m_Image = mitk::IOUtil::LoadImage(imagePath); m_ReferenceImage = NULL; - m_DenoisedImage = mitk::DiffusionImage::New(); + m_DenoisedImage = mitk::Image::New(); //initialise Filter m_DenoisingFilter = itk::NonLocalMeansDenoisingFilter::New(); - m_DenoisingFilter->SetInputImage(m_Image->GetVectorImage()); -// m_DenoisingFilter->SetInputMask(m_ImageMask); + VectorImagetType::Pointer vectorImage; + mitk::CastToItkImage(m_Image,vectorImage); + m_DenoisingFilter->SetInputImage(vectorImage); m_DenoisingFilter->SetNumberOfThreads(1); m_DenoisingFilter->SetComparisonRadius(1); m_DenoisingFilter->SetSearchRadius(1); m_DenoisingFilter->SetVariance(500); } void tearDown() { m_Image = NULL; m_ImageMask = NULL; m_ReferenceImage = NULL; m_DenoisingFilter = NULL; m_DenoisedImage = NULL; } void Denoise_NLMg_shouldReturnTrue() { std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMg.dwi"); - m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); + m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath); m_DenoisingFilter->SetUseRicianAdaption(false); m_DenoisingFilter->SetUseJointInformation(false); try { m_DenoisingFilter->Update(); } catch(std::exception& e) { MITK_ERROR << e.what(); } - m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); - m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); - m_DenoisedImage->SetDirections(m_Image->GetDirections()); - m_DenoisedImage->InitializeFromVectorImage(); + mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage); + m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone()); MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMg should always return the same result."); } void Denoise_NLMr_shouldReturnTrue() { std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMr.dwi"); - m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); + m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath); m_DenoisingFilter->SetUseRicianAdaption(true); m_DenoisingFilter->SetUseJointInformation(false); try { m_DenoisingFilter->Update(); } catch(std::exception& e) { MITK_ERROR << e.what(); } - m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); - m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); - m_DenoisedImage->SetDirections(m_Image->GetDirections()); - m_DenoisedImage->InitializeFromVectorImage(); + mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage); + m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone()); MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMr should always return the same result."); } void Denoise_NLMv_shouldReturnTrue() { std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMv.dwi"); - m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); - + m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath); m_DenoisingFilter->SetUseRicianAdaption(false); m_DenoisingFilter->SetUseJointInformation(true); try { m_DenoisingFilter->Update(); } catch(std::exception& e) { MITK_ERROR << e.what(); } - m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); - m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); - m_DenoisedImage->SetDirections(m_Image->GetDirections()); - m_DenoisedImage->InitializeFromVectorImage(); + mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage); + m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone()); MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMv should always return the same result."); } void Denoise_NLMvr_shouldReturnTrue() { std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMvr.dwi"); - m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); + m_ReferenceImage = mitk::IOUtil::LoadImage(referenceImagePath); m_DenoisingFilter->SetUseRicianAdaption(true); m_DenoisingFilter->SetUseJointInformation(true); try { m_DenoisingFilter->Update(); } catch(std::exception& e) { MITK_ERROR << e.what(); } - m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); - m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); - m_DenoisedImage->SetDirections(m_Image->GetDirections()); - m_DenoisedImage->InitializeFromVectorImage(); + mitk::GrabItkImageMemory(m_DenoisingFilter->GetOutput(),m_DenoisedImage); + m_DenoisedImage->SetPropertyList(m_Image->GetPropertyList()->Clone()); MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMvr should always return the same result."); } }; MITK_TEST_SUITE_REGISTRATION(mitkNonLocalMeansDenoising) diff --git a/Modules/DiffusionImaging/DiffusionCore/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake index d5ad2057eb..dae157948b 100644 --- a/Modules/DiffusionImaging/DiffusionCore/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake @@ -1,150 +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 DicomImport/mitkDiffusionDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderSiemensDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderSiemensDICOMFileHelper.cpp DicomImport/mitkDiffusionHeaderSiemensMosaicDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderGEDICOMFileReader.cpp DicomImport/mitkDiffusionHeaderPhilipsDICOMFileReader.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp - IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp - - IODataStructures/DiffusionWeightedImages/mitkImageToDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCorrectionFilter.cpp # Properties IODataStructures/Properties/mitkBValueMapProperty.cpp IODataStructures/Properties/mitkGradientDirectionsProperty.cpp IODataStructures/Properties/mitkMeasurementFrameProperty.cpp IODataStructures/Properties/mitkDiffusionPropertyHelper.cpp IODataStructures/Properties/mitkNodePredicateIsDWI.cpp # Serializer IODataStructures/Properties/mitkBValueMapPropertySerializer.cpp IODataStructures/Properties/mitkGradientDirectionsPropertySerializer.cpp IODataStructures/Properties/mitkMeasurementFramePropertySerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkQBallImage.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImage.cpp Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/itkDwiGradientLengthCorrectionFilter.cpp Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h # Registration Algorithms & Co. Algorithms/Registration/mitkRegistrationWrapper.cpp Algorithms/Registration/mitkPyramidImageRegistrationMethod.cpp # Algorithms/Registration/mitkRegistrationMethodITK4.cpp + Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.cpp # MultishellProcessing Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.cpp Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h # Rendering - Rendering/mitkDiffusionImageMapper.h Rendering/mitkOdfVtkMapper2D.h # Reconstruction Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Algorithms/Reconstruction/itkPointShell.h Algorithms/Reconstruction/itkOrientationDistributionFunction.h Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h # MultishellProcessing Algorithms/Reconstruction/MultishellProcessing/itkRadialMultishellToSingleshellImageFilter.h Algorithms/Reconstruction/MultishellProcessing/itkDWIVoxelFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.h Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.h - # IO Datastructures - IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h + # Properties + IODataStructures/Properties/mitkBValueMapProperty.h + IODataStructures/Properties/mitkGradientDirectionsProperty.h + IODataStructures/Properties/mitkMeasurementFrameProperty.h + IODataStructures/Properties/mitkDiffusionPropertyHelper.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h Algorithms/itkMergeDiffusionImagesFilter.h Algorithms/itkDwiPhantomGenerationFilter.h Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h Algorithms/itkMrtrixPeakImageConverter.h Algorithms/itkFslPeakImageConverter.h Algorithms/itkShCoefficientImageImporter.h Algorithms/itkShCoefficientImageExporter.h Algorithms/itkOdfMaximaExtractionFilter.h Algorithms/itkResampleDwiImageFilter.h Algorithms/itkDwiGradientLengthCorrectionFilter.h Algorithms/itkAdcImageFilter.h Algorithms/itkDwiNormilzationFilter.h Algorithms/itkSplitDWImageFilter.h Algorithms/itkRemoveDwiChannelFilter.h Algorithms/itkExtractDwiChannelFilter.h Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h - Algorithms/mitkDiffusionImageToDiffusionImageFilter.h Algorithms/itkNonLocalMeansDenoisingFilter.h Algorithms/itkVectorImageToImageFilter.h ) set( TOOL_FILES )