diff --git a/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt b/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt index 9fc7c8ab32..602732e5a1 100644 --- a/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt +++ b/Modules/DiffusionImaging/DiffusionCore/CMakeLists.txt @@ -1,21 +1,21 @@ # With apple gcc 4.2.1 the following waring leads to an build error if boost is enabled if(APPLE) mitkFunctionCheckCAndCXXCompilerFlags("-Wno-error=empty-body" CMAKE_C_FLAGS CMAKE_CXX_FLAGS) endif() MITK_CREATE_MODULE( SUBPROJECTS MITK-DTI - INCLUDE_DIRS include/ include/Algorithms include/Algorithms/Reconstruction include/Algorithms/Registration include/Algorithms/Reconstruction/MultishellProcessing include/DicomImport include/IODataStructures/DiffusionWeightedImages include/IODataStructures/Properties include/IODataStructures/QBallImages include/IODataStructures/TensorImages include/IODataStructures include/Rendering ${CMAKE_CURRENT_BINARY_DIR} + INCLUDE_DIRS include/ include/Algorithms include/Algorithms/Reconstruction include/Algorithms/Registration include/Algorithms/Reconstruction/MultishellProcessing include/Algorithms/Reconstruction/FittingFunctions include/DicomImport include/IODataStructures/DiffusionWeightedImages include/IODataStructures/Properties include/IODataStructures/QBallImages include/IODataStructures/TensorImages include/IODataStructures include/Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS MitkMapperExt MitkPlanarFigure MitkImageExtraction MitkSceneSerializationBase MitkDICOMReader PACKAGE_DEPENDS PUBLIC ITK|ITKTestKernel+ITKRegistrationCommon+ITKMetricsv4+ITKRegistrationMethodsv4+ITKDistanceMap+ITKLabelVoting+ITKVTK PUBLIC VTK|vtkFiltersProgrammable ) if(MSVC) mitkFunctionCheckCAndCXXCompilerFlags("/wd4005" CMAKE_C_FLAGS CMAKE_CXX_FLAGS) endif() add_subdirectory(Testing) add_subdirectory(cmdapps) add_subdirectory(autoload/IO) diff --git a/Modules/DiffusionImaging/DiffusionCore/files.cmake b/Modules/DiffusionImaging/DiffusionCore/files.cmake index 7c400740d8..c20460b998 100644 --- a/Modules/DiffusionImaging/DiffusionCore/files.cmake +++ b/Modules/DiffusionImaging/DiffusionCore/files.cmake @@ -1,150 +1,157 @@ set(CPP_FILES # DicomImport # 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/mitkDiffusionImageCorrectionFilter.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageCreationFilter.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 # 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 include/mitkDiffusionFunctionCollection.h # Rendering include/Rendering/mitkOdfVtkMapper2D.h # Reconstruction include/Algorithms/Reconstruction/itkDiffusionQballReconstructionImageFilter.h include/Algorithms/Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h include/Algorithms/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h include/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h include/Algorithms/Reconstruction/itkPointShell.h include/Algorithms/Reconstruction/itkOrientationDistributionFunction.h include/Algorithms/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h include/Algorithms/Reconstruction/itkDiffusionKurtosisReconstructionImageFilter.h + include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.h + include/Algorithms/Reconstruction/itkMultiTensorImageFilter.h + + # Fitting functions + include/Algorithms/Reconstruction/FittingFunctions/mitkAbstractFitter.h + include/Algorithms/Reconstruction/FittingFunctions/mitkMultiTensorFitter.h + include/Algorithms/Reconstruction/FittingFunctions/mitkBallStickFitter.h + # MultishellProcessing include/Algorithms/Reconstruction/MultishellProcessing/itkRadialMultishellToSingleshellImageFilter.h include/Algorithms/Reconstruction/MultishellProcessing/itkDWIVoxelFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkADCAverageFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkKurtosisFitFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkBiExpFitFunctor.h include/Algorithms/Reconstruction/MultishellProcessing/itkADCFitFunctor.h # Properties include/IODataStructures/Properties/mitkBValueMapProperty.h include/IODataStructures/Properties/mitkGradientDirectionsProperty.h include/IODataStructures/Properties/mitkMeasurementFrameProperty.h include/IODataStructures/Properties/mitkDiffusionPropertyHelper.h include/IODataStructures/DiffusionWeightedImages/mitkDiffusionImageTransformedCreationFilter.h # Algorithms include/Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h include/Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h include/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h include/Algorithms/itkTensorDerivedMeasurementsFilter.h include/Algorithms/itkBrainMaskExtractionImageFilter.h include/Algorithms/itkB0ImageExtractionImageFilter.h include/Algorithms/itkB0ImageExtractionToSeparateImageFilter.h include/Algorithms/itkTensorImageToDiffusionImageFilter.h include/Algorithms/itkTensorToL2NormImageFilter.h include/Algorithms/itkGaussianInterpolateImageFunction.h include/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h include/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h include/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.h include/Algorithms/itkCartesianToPolarVectorImageFilter.h include/Algorithms/itkPolarToCartesianVectorImageFilter.h include/Algorithms/itkDistanceMapFilter.h include/Algorithms/itkProjectionFilter.h include/Algorithms/itkResidualImageFilter.h include/Algorithms/itkExtractChannelFromRgbaImageFilter.h include/Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h include/Algorithms/itkMergeDiffusionImagesFilter.h - include/Algorithms/itkDwiPhantomGenerationFilter.h include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.h include/Algorithms/itkShCoefficientImageImporter.h include/Algorithms/itkShCoefficientImageExporter.h include/Algorithms/itkOdfMaximaExtractionFilter.h include/Algorithms/itkResampleDwiImageFilter.h include/Algorithms/itkDwiGradientLengthCorrectionFilter.h include/Algorithms/itkAdcImageFilter.h include/Algorithms/itkDwiNormilzationFilter.h include/Algorithms/itkSplitDWImageFilter.h include/Algorithms/itkRemoveDwiChannelFilter.h include/Algorithms/itkExtractDwiChannelFilter.h include/Algorithms/Registration/mitkDWIHeadMotionCorrectionFilter.h include/Algorithms/itkNonLocalMeansDenoisingFilter.h include/Algorithms/itkVectorImageToImageFilter.h ) set( TOOL_FILES ) diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkAbstractFitter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkAbstractFitter.h new file mode 100644 index 0000000000..bed905a605 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkAbstractFitter.h @@ -0,0 +1,124 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef _MITK_AbstractFitter_H +#define _MITK_AbstractFitter_H + +#include +#include +#include + +namespace mitk { + +class AbstractFitter: public vnl_least_squares_function +{ + +public: + typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientContainerType; + typedef itk::VectorImage< unsigned short, 3 > DiffusionImageType; + typedef itk::DiffusionTensor3D TensorType; + + AbstractFitter(unsigned int number_of_parameters, unsigned int number_of_measurements) : + vnl_least_squares_function(number_of_parameters, number_of_measurements, no_gradient) + { + + } + + GradientContainerType::Pointer gradientDirections; + typename DiffusionImageType::PixelType measurements; + std::vector bValues; + double S0; + std::vector weightedIndices; + + void set_S0(double val) + { + S0 = val; + } + + void set_measurements(const typename DiffusionImageType::PixelType& m) + { + measurements = m; + } + + void set_bvalues(std::vector& x) + { + bValues = x; + } + + void set_weightedIndices(std::vector& x) + { + weightedIndices = x; + } + + void set_gradient_directions(const GradientContainerType::Pointer& directions) + { + gradientDirections = directions; + } + + static void Sph2Cart(vnl_vector_fixed& dir, const double &theta, const double &phi) + { + dir[0] = std::sin(theta)*std::cos(phi); + dir[1] = std::sin(theta)*std::sin(phi); + dir[2] = std::cos(theta); + } + + static void Cart2Sph(const vnl_vector_fixed& dir, double &theta, double &phi) + { + theta = std::acos( dir[2] ); + + // phi goes from 0.0 (+x axis) and wraps at 2 * PI + // theta goes from 0.0 (+z axis) and wraps at PI + + // if x and y are 0.0 or very close, return phi == 0 + if (fabs(dir[0]) + fabs(dir[1]) < 1E-9) + { + phi = 0.0; + } + else { + // ie, if ( x == 0 && y == 0 ) == false + if (dir[1] == 0.0) + { + if (dir[0] > 0.0) + phi = 0.0; + else + phi = M_PI; + } + else if (dir[0] == 0.0) + { + // avoid div by zero + if (dir[1] > 0) { + phi = M_PI / 2.0; + } + else { + phi = 1.5 * M_PI; + } + } + else if (dir[0] > 0.0 && dir[1] > 0.0) // first quadrant + phi = atan(dir[1] / dir[0]); + else if (dir[0] < 0.0 && dir[1] > 0.0) // second quadrant + phi = M_PI + atan(dir[1] / dir[0]); + else if (dir[0] < 0.0 && dir[1] < 0.0) // third quadrant + phi = M_PI + atan(dir[1] / dir[0]); + else // fourth quadrant + phi = 2.0 * M_PI + atan(dir[1] / dir[0]); + } + } + +}; +} + +#endif diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkBallStickFitter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkBallStickFitter.h new file mode 100644 index 0000000000..7e9acd51a6 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkBallStickFitter.h @@ -0,0 +1,72 @@ +#include +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef _MITK_BallStickFitter_H +#define _MITK_BallStickFitter_H + +#include + +namespace mitk { + +struct BallStickFitter: public AbstractFitter +{ +public : + + BallStickFitter(unsigned int number_of_parameters, unsigned int number_of_measurements) : + AbstractFitter(number_of_parameters, number_of_measurements) + { + + } + + double penalty(const vnl_vector& x) + { + double p = 0; + + if (x[0]<0 || x[0]>1) + p += 10e6; + if (x[1]<0) + p += 10e6; + + return p; + } + + void f(const vnl_vector& x, vnl_vector& fx) override + { + const double & f = x[0]; + const double & d = x[1]; + const double & theta = x[2]; + const double & phi = x[3]; + vnl_vector_fixed dir; + Sph2Cart(dir, theta, phi); + + for(auto s : weightedIndices) + { + double s_iso = S0 * std::exp(-bValues[s] * d); + GradientDirectionType g = gradientDirections->GetElement(s); + g.normalize(); + double dot = dot_product(g, dir); + double s_aniso = S0 * std::exp(-bValues[s] * d * dot*dot ); + double approx = (1-f)*s_iso + f*s_aniso; + const double factor = measurements[s] - approx; + fx[s] = factor*factor + penalty(x); + } + } +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkMultiTensorFitter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkMultiTensorFitter.h new file mode 100644 index 0000000000..7018e02f20 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/FittingFunctions/mitkMultiTensorFitter.h @@ -0,0 +1,117 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef _MITK_MultiTensorFitter_H +#define _MITK_MultiTensorFitter_H + +#include + +namespace mitk { + +class MultiTensorFitter: public AbstractFitter +{ +public: + + MultiTensorFitter(unsigned int number_of_tensors=2, unsigned int number_of_measurements=1) : + AbstractFitter(check(number_of_tensors), number_of_measurements) + { + num_tensors = number_of_tensors; + } + + int num_tensors; + + int check(int n) + { + if (n==1) + return 6; + return n*7; + } + + vnl_vector_fixed GetLargestEv(TensorType& tensor) + { + TensorType::EigenValuesArrayType eigenvalues; + TensorType::EigenVectorsMatrixType eigenvectors; + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + vnl_vector_fixed ev; + ev[0] = eigenvectors(2, 0); + ev[1] = eigenvectors(2, 1); + ev[2] = eigenvectors(2, 2); + if (ev.magnitude()>mitk::eps) + ev.normalize(); + else + ev.fill(0.0); + + return ev; + } + + double penalty(const vnl_vector& x) + { + double p = 0; + + if (num_tensors>1) + { + double f = 0; + for (int i=0; i& x, vnl_vector& fx) override + { + int elements = 7; + for(auto s : weightedIndices) + { + GradientDirectionType g = gradientDirections->GetElement(s); + g.normalize(); + + itk::DiffusionTensor3D< double > S; + S[0] = g[0]*g[0]; + S[1] = g[1]*g[0]; + S[2] = g[2]*g[0]; + S[3] = g[1]*g[1]; + S[4] = g[2]*g[1]; + S[5] = g[2]*g[2]; + + double approx = 0; + for (int i=0; i1 && x.size()>6) + signal *= x[elements-1+i*elements]; + approx += signal; + } + + const double factor = measurements[s] - approx; + + fx[s] = factor*factor + penalty(x); + } + } + +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.h new file mode 100644 index 0000000000..d8f86031e5 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.h @@ -0,0 +1,103 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef __itkBallAndSticksImageFilter_h_ +#define __itkBallAndSticksImageFilter_h_ + +#include "itkImageToImageFilter.h" +#include "itkVectorImage.h" +#include +#include +#include +#include +#include +#include + +#define NUM_TENSORS 2 + +namespace itk{ +/** \class BallAndSticksImageFilter + */ + +template< class TInPixelType, class TOutPixelType > +class BallAndSticksImageFilter : + public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > +{ + +public: + + typedef BallAndSticksImageFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > Superclass; + typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientContainerType; + typedef itk::Image< unsigned char, 3> ItkUcharImageType; + typedef itk::Image< float, 4 > PeakImageType; + + typedef itk::DiffusionTensor3DReconstructionImageFilter< TInPixelType, TInPixelType, float > TensorRecFilterType; + typedef itk::DiffusionTensor3D TensorType; + typedef itk::Image TensorImageType; + + /** Method for creation through the object factory. */ + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(BallAndSticksImageFilter, ImageToImageFilter) + + typedef typename Superclass::InputImageType InputImageType; + typedef typename Superclass::OutputImageType OutputImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + + itkSetMacro( MaskImage, ItkUcharImageType::Pointer ) + itkSetMacro( B_value, double ) + itkSetMacro( GradientDirections, GradientContainerType::Pointer ) + itkGetMacro( PeakImage, PeakImageType::Pointer ) + itkGetMacro( OutDwi, typename InputImageType::Pointer ) + + protected: + BallAndSticksImageFilter(); + ~BallAndSticksImageFilter() {} + void PrintSelf(std::ostream& os, Indent indent) const; + + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); + + + double m_B_value; + std::vector m_B_values; + std::vector m_WeightedIndices; + std::vector m_UnWeightedIndices; + GradientContainerType::Pointer m_GradientDirections; + + ItkUcharImageType::Pointer m_MaskImage; + PeakImageType::Pointer m_PeakImage; + TensorImageType::Pointer m_TensorImage; + typename InputImageType::Pointer m_OutDwi; + + vnl_vector FitSingleVoxel( const typename InputImageType::PixelType &input, const typename InputImageType::IndexType &idx); + +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkBallAndSticksImageFilter.txx" +#endif + +#endif //__itkBallAndSticksImageFilter_h_ + diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.txx new file mode 100644 index 0000000000..c708ae2170 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkBallAndSticksImageFilter.txx @@ -0,0 +1,252 @@ +/*=================================================================== + +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 __itkBallAndSticksImageFilter_txx +#define __itkBallAndSticksImageFilter_txx + +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionIterator.h" +#include + + +namespace itk { + + +template< class TInPixelType, class TOutPixelType > +BallAndSticksImageFilter< TInPixelType, TOutPixelType> +::BallAndSticksImageFilter() + : m_B_value(0) +{ + this->SetNumberOfRequiredInputs( 1 ); +} + +template< class TInPixelType, class TOutPixelType > +void +BallAndSticksImageFilter< TInPixelType, TOutPixelType> +::BeforeThreadedGenerateData() +{ + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + outputImage->FillBuffer(0.0); + + m_UnWeightedIndices.clear(); + m_WeightedIndices.clear(); + m_B_values.clear(); + + for (unsigned int i=0; iSize(); i++) + { + GradientDirectionType g = m_GradientDirections->GetElement(i); + double twonorm = g.two_norm(); + double b = m_B_value*twonorm*twonorm; + m_B_values.push_back(b); + + if (b>100) + m_WeightedIndices.push_back(i); + else + m_UnWeightedIndices.push_back(i); + } + + if (m_UnWeightedIndices.empty()) + mitkThrow() << "Unweighted (b=0 s/mm²) image volume missing!"; + + + itk::Vector spacing3 = outputImage->GetSpacing(); + itk::Point origin3 = outputImage->GetOrigin(); + itk::Matrix direction3 = outputImage->GetDirection(); + itk::ImageRegion<3> imageRegion3 = outputImage->GetLargestPossibleRegion(); + + itk::Vector spacing4; + itk::Point origin4; + itk::Matrix direction4; + itk::ImageRegion<4> imageRegion4; + + spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; + origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin4[3] = 0; + for (int r=0; r<3; r++) + for (int c=0; c<3; c++) + direction4[r][c] = direction3[r][c]; + direction4[3][3] = 1; + imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); + imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); + imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); + imageRegion4.SetSize(3, 3); + + m_PeakImage = PeakImageType::New(); + m_PeakImage->SetSpacing( spacing4 ); + m_PeakImage->SetOrigin( origin4 ); + m_PeakImage->SetDirection( direction4 ); + m_PeakImage->SetRegions( imageRegion4 ); + m_PeakImage->Allocate(); + m_PeakImage->FillBuffer(0.0); + + m_OutDwi = InputImageType::New(); + m_OutDwi->SetSpacing( spacing3 ); + m_OutDwi->SetOrigin( origin3 ); + m_OutDwi->SetDirection( direction3 ); + m_OutDwi->SetRegions( imageRegion3 ); + m_OutDwi->SetVectorLength(m_GradientDirections->Size()); + m_OutDwi->Allocate(); + + MITK_INFO << "Initial tensor fit"; + typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); + typename TensorRecFilterType::Pointer tensorReconstructionFilter = TensorRecFilterType::New(); + tensorReconstructionFilter->SetBValue( m_B_value ); + tensorReconstructionFilter->SetGradientImage( m_GradientDirections, inputImagePointer ); + tensorReconstructionFilter->Update(); + m_TensorImage = tensorReconstructionFilter->GetOutput(); +} + +template< class TInPixelType, class TOutPixelType > +vnl_vector +BallAndSticksImageFilter< TInPixelType, TOutPixelType>::FitSingleVoxel( const typename InputImageType::PixelType &input, const typename InputImageType::IndexType &idx) +{ + double S0 = 0; + for (auto i : m_UnWeightedIndices) + S0 += input[i]; + S0 /= m_UnWeightedIndices.size(); + + // Linear tensor fit + double md = 0.001; + double f = 0.5; + double theta = 0; + double phi = 0; + { + TensorType tensor = m_TensorImage->GetPixel(idx); + TensorType::EigenValuesArrayType eigenvalues; + TensorType::EigenVectorsMatrixType eigenvectors; + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + vnl_vector_fixed ev; + ev[0] = eigenvectors(2, 0); + ev[1] = eigenvectors(2, 1); + ev[2] = eigenvectors(2, 2); + if (ev.magnitude()>mitk::eps) + ev.normalize(); + else + ev.fill(0.0); + + md = fabs(eigenvalues[0]+eigenvalues[1]+eigenvalues[2])/3; + f = tensor.GetFractionalAnisotropy(); + if (f<0.1) + f = 0.1; + else if (f>0.9) + f = 0.9; + + mitk::AbstractFitter::Cart2Sph(ev, theta, phi); + } + + int num_params = 4; // f, d, phi, theta + + vnl_vector x; x.set_size(num_params); x.fill(0.0); + x[0] = f; // f (volume fraction) + x[1] = md; // d + x[2] = theta; + x[3] = phi; + + mitk::BallStickFitter bs_fit(num_params, m_GradientDirections->size()); + bs_fit.set_S0(S0); + bs_fit.set_weightedIndices(m_WeightedIndices); + bs_fit.set_bvalues(m_B_values); + bs_fit.set_gradient_directions(m_GradientDirections); + bs_fit.set_measurements(input); + + vnl_levenberg_marquardt lm(bs_fit); + lm.minimize(x); + return x; +} + +template< class TInPixelType, class TOutPixelType > +void +BallAndSticksImageFilter< TInPixelType, TOutPixelType> +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType ) +{ + typename OutputImageType::Pointer outputImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + + ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); + oit.GoToBegin(); + + typedef ImageRegionConstIterator< InputImageType > InputIteratorType; + typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); + + InputIteratorType git( inputImagePointer, outputRegionForThread ); + git.GoToBegin(); + while( !git.IsAtEnd() ) + { + typename InputImageType::PixelType pix = git.Get(); + vnl_vector x = FitSingleVoxel(pix, git.GetIndex()); + + /// PEAKS + PeakImageType::IndexType idx4; + idx4[0] = oit.GetIndex()[0]; + idx4[1] = oit.GetIndex()[1]; + idx4[2] = oit.GetIndex()[2]; + + oit.Set( x[0] ); + vnl_vector_fixed dir; + mitk::AbstractFitter::Sph2Cart(dir, x[2], x[3]); + + idx4[3] = 0; + m_PeakImage->SetPixel(idx4, dir[0]); + idx4[3] = 1; + m_PeakImage->SetPixel(idx4, dir[1]); + idx4[3] = 2; + m_PeakImage->SetPixel(idx4, dir[2]); + + /// DWI from ball-stick + typename InputImageType::PixelType dPix; dPix.SetSize(m_GradientDirections->Size()); dPix.Fill(0.0); + for (unsigned int i=0; iSize(); i++) + { + GradientDirectionType g = m_GradientDirections->GetElement(i); + double twonorm = g.two_norm(); + double b = m_B_value*twonorm*twonorm; + g.normalize(); + + double s_iso = 1000 * std::exp(-b * x[1]); + + double dot = dot_product(g, dir); + double s_aniso = 1000 * std::exp(-b * x[1] * dot*dot ); + + double approx = (1-x[0])*s_iso + x[0]*s_aniso; + + dPix[i] = approx; + } + m_OutDwi->SetPixel(oit.GetIndex(), dPix); + + ++oit; + ++git; + } + + std::cout << "One Thread finished calculation" << std::endl; +} + +template< class TInPixelType, class TOutPixelType > +void +BallAndSticksImageFilter< TInPixelType, TOutPixelType> +::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os,indent); +} + +} + +#endif // __itkBallAndSticksImageFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkMultiTensorImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkMultiTensorImageFilter.h new file mode 100644 index 0000000000..4dff86d146 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkMultiTensorImageFilter.h @@ -0,0 +1,198 @@ +/*=================================================================== + +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 __itkMultiTensorImageFilter_h_ +#define __itkMultiTensorImageFilter_h_ + +#include "itkImageToImageFilter.h" +#include "itkVectorImage.h" +#include +#include +#include + +namespace itk{ +/** \class MultiTensorImageFilter + */ + +template< class TInPixelType, class TOutPixelType > +class MultiTensorImageFilter : + public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > +{ + +public: + + typedef MultiTensorImageFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > Superclass; + typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientContainerType; + typedef itk::Image< unsigned char, 3> ItkUcharImageType; + typedef itk::Image< float, 4 > PeakImageType; + typedef itk::DiffusionTensor3D TensorType; + typedef itk::Image TensorImageType; + + /** Method for creation through the object factory. */ + itkFactorylessNewMacro(Self) + itkCloneMacro(Self) + + /** Runtime information support. */ + itkTypeMacro(MultiTensorImageFilter, ImageToImageFilter) + + typedef typename Superclass::InputImageType InputImageType; + typedef typename Superclass::OutputImageType OutputImageType; + typedef typename Superclass::OutputImageRegionType OutputImageRegionType; + + itkSetMacro( MaskImage, ItkUcharImageType::Pointer ) + itkSetMacro( B_value, double ) + itkSetMacro( GradientDirections, GradientContainerType::Pointer ) + + itkGetMacro( PeakImage, PeakImageType::Pointer ) + std::vector< TensorImageType::Pointer > GetTensorImages(){ return m_TensorImages; } + +protected: + MultiTensorImageFilter(); + ~MultiTensorImageFilter() {} + void PrintSelf(std::ostream& os, Indent indent) const; + + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); + + double m_B_value; + std::vector m_B_values; + std::vector m_WeightedIndices; + std::vector m_UnWeightedIndices; + + GradientContainerType::Pointer m_GradientDirections; + + ItkUcharImageType::Pointer m_MaskImage; + PeakImageType::Pointer m_PeakImage; + + std::vector< TensorImageType::Pointer > m_TensorImages; + int m_NumTensors; + + vnl_vector FitSingleVoxel( const typename InputImageType::PixelType &input); + +// struct multiTensorLeastSquaresFunction: public vnl_least_squares_function +// { +// GradientContainerType::Pointer gradientDirections; +// typename InputImageType::PixelType measurements; +// std::vector bValues; +// double S0; +// int num_tensors; +// std::vector weightedIndices; + +// void set_S0(double val) +// { +// S0 = val; +// } + +// void set_measurements(const typename InputImageType::PixelType& m) +// { +// measurements = m; +// } + +// void set_bvalues(std::vector& x) +// { +// bValues = x; +// } + +// void set_weightedIndices(std::vector& x) +// { +// weightedIndices = x; +// } + +// void set_gradient_directions(const GradientContainerType::Pointer& directions) +// { +// gradientDirections = directions; +// } + +// int check(int n) +// { +// if (n==1) +// return 6; +// return n*7; +// } + +// multiTensorLeastSquaresFunction(unsigned int number_of_tensors=2, unsigned int number_of_measurements=1) : +// vnl_least_squares_function(check(number_of_tensors), number_of_measurements, no_gradient) +// { +// num_tensors = number_of_tensors; +// } + +// void f(const vnl_vector& x, vnl_vector& fx) override +// { +// int elements = 7; +// for(auto s : weightedIndices) +// { +// GradientDirectionType g = gradientDirections->GetElement(s); +// g.normalize(); + +// itk::DiffusionTensor3D< double > S; +// S[0] = g[0]*g[0]; +// S[1] = g[1]*g[0]; +// S[2] = g[2]*g[0]; +// S[3] = g[1]*g[1]; +// S[4] = g[2]*g[1]; +// S[5] = g[2]*g[2]; + +// double approx = 0; +// double penalty = 0; +// double f = 0; +// if (num_tensors==1 && x.size()==6) +// f = 1; + +// for (int i=0; i1 && x.size()>6 && x[elements-1+i*elements]<0) // penalty if volume fractiuon is < 0 +// penalty += 10e6; + +// if (num_tensors>1 && x.size()>6) +// f += fabs(x[elements-1+i*elements]); + +// double D = x[0+i*elements]*S[0] + x[1+i*elements]*S[1] + x[2+i*elements]*S[2] + +// x[1+i*elements]*S[1] + x[3+i*elements]*S[3] + x[4+i*elements]*S[4] + +// x[2+i*elements]*S[2] + x[4+i*elements]*S[4] + x[5+i*elements]*S[5]; +// double signal = S0 * std::exp ( -bValues[s] * D ); +// if (num_tensors>1 && x.size()>6) +// signal *= x[elements-1+i*elements]; +// approx += signal; +// } + +// const double factor = measurements[s] - approx; +// penalty += fabs(f-1)*10e7; // penalty if volume fractiuon sum is != 1 + +// fx[s] = factor*factor + penalty; +// } +// } + +// }; + +}; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkMultiTensorImageFilter.txx" +#endif + +#endif //__itkMultiTensorImageFilter_h_ + diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkMultiTensorImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkMultiTensorImageFilter.txx new file mode 100644 index 0000000000..18430b4405 --- /dev/null +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/Reconstruction/itkMultiTensorImageFilter.txx @@ -0,0 +1,253 @@ +/*=================================================================== + +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 __itkMultiTensorImageFilter_txx +#define __itkMultiTensorImageFilter_txx + +#include +#include +#include + +#define _USE_MATH_DEFINES +#include + +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionIterator.h" +#include +#include + +namespace itk { + + +template< class TInPixelType, class TOutPixelType > +MultiTensorImageFilter< TInPixelType, TOutPixelType> +::MultiTensorImageFilter() + : m_B_value(0) + , m_NumTensors(2) +{ + this->SetNumberOfRequiredInputs( 1 ); +} + +template< class TInPixelType, class TOutPixelType > +void +MultiTensorImageFilter< TInPixelType, TOutPixelType> +::BeforeThreadedGenerateData() +{ + typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + outputImage->FillBuffer(0.0); + + m_UnWeightedIndices.clear(); + m_WeightedIndices.clear(); + m_B_values.clear(); + + for (unsigned int i=0; iSize(); i++) + { + GradientDirectionType g = m_GradientDirections->GetElement(i); + double twonorm = g.two_norm(); + double b = m_B_value*twonorm*twonorm; + m_B_values.push_back(b); + + if (b>100) + m_WeightedIndices.push_back(i); + else + m_UnWeightedIndices.push_back(i); + } + + if (m_UnWeightedIndices.empty()) + mitkThrow() << "Unweighted (b=0 s/mm²) image volume missing!"; + + itk::Vector spacing3 = outputImage->GetSpacing(); + itk::Point origin3 = outputImage->GetOrigin(); + itk::Matrix direction3 = outputImage->GetDirection(); + itk::ImageRegion<3> imageRegion3 = outputImage->GetLargestPossibleRegion(); + + itk::Vector spacing4; + itk::Point origin4; + itk::Matrix direction4; + itk::ImageRegion<4> imageRegion4; + + spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; + origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin4[3] = 0; + for (int r=0; r<3; r++) + for (int c=0; c<3; c++) + direction4[r][c] = direction3[r][c]; + direction4[3][3] = 1; + imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); + imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); + imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); + imageRegion4.SetSize(3, 3*m_NumTensors); + + m_PeakImage = PeakImageType::New(); + m_PeakImage->SetSpacing( spacing4 ); + m_PeakImage->SetOrigin( origin4 ); + m_PeakImage->SetDirection( direction4 ); + m_PeakImage->SetRegions( imageRegion4 ); + m_PeakImage->Allocate(); + m_PeakImage->FillBuffer(0.0); + + m_TensorImages.clear(); + for (int i=0; iSetSpacing( spacing3 ); + tImg->SetOrigin( origin3 ); + tImg->SetDirection( direction3 ); + tImg->SetRegions( imageRegion3 ); + tImg->Allocate(); + m_TensorImages.push_back(tImg); + } +} + +template< class TInPixelType, class TOutPixelType > +vnl_vector +MultiTensorImageFilter< TInPixelType, TOutPixelType>::FitSingleVoxel( const typename InputImageType::PixelType &input) +{ + double S0 = 0; + for (auto i : m_UnWeightedIndices) + S0 += input[i]; + S0 /= m_UnWeightedIndices.size(); + + // signle tensor fit + mitk::MultiTensorFitter ls_fit(1, m_GradientDirections->size()); + ls_fit.set_S0(S0); + ls_fit.set_weightedIndices(m_WeightedIndices); + ls_fit.set_bvalues(m_B_values); + ls_fit.set_gradient_directions(m_GradientDirections); + ls_fit.set_measurements(input); + + vnl_levenberg_marquardt lm(ls_fit); + vnl_vector x; + x.set_size(6); + x.fill(0.0); + + lm.minimize(x); + +// TensorType tensor; + + vnl_vector y; + y.set_size(7*m_NumTensors); + y.fill(0.0); + for (int i=0; i<6; i++) + { + y[i] = x[i]; +// tensor[i] = x[i]; + } + y[6]=1; +// double fa = tensor.GetFractionalAnisotropy(); + + if (m_NumTensors>1) + { + mitk::MultiTensorFitter ls_fit(m_NumTensors, m_GradientDirections->size()); + ls_fit.set_S0(S0); + ls_fit.set_weightedIndices(m_WeightedIndices); + ls_fit.set_bvalues(m_B_values); + ls_fit.set_gradient_directions(m_GradientDirections); + ls_fit.set_measurements(input); + + vnl_levenberg_marquardt lm(ls_fit); + + for (int i=0; i +void +MultiTensorImageFilter< TInPixelType, TOutPixelType> +::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType ) +{ + typedef ImageRegionConstIterator< InputImageType > InputIteratorType; + typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); + + InputIteratorType git( inputImagePointer, outputRegionForThread ); + git.GoToBegin(); + while( !git.IsAtEnd() ) + { + typename InputImageType::PixelType pix = git.Get(); + vnl_vector x = FitSingleVoxel(pix); + + typedef itk::DiffusionTensor3D TensorType; + + int elements = 7; + for (int t=0; tSetPixel(git.GetIndex(), tensor); + } + +// TensorType tensor; +// tensor.Fill(0.0); +// tensor[0] = x[0]*x[1]; +// tensor[3] = tensor[0]; +// tensor[5] = tensor[0]; +// m_TensorImages.at(m_NumTensors+1)->SetPixel(git.GetIndex(), tensor); + +// /// DWI from tensor +// int elements = 7; +// int num_ten = m_NumTensors; +// typename DiffusionImageType::PixelType dPix; dPix.SetSize(m_GradientDirections->Size()); dPix.Fill(0.0); +// for (unsigned int i=0; iSize(); i++) +// { +// GradientDirectionType g = m_GradientDirections->GetElement(i); +// double twonorm = g.two_norm(); +// double b = m_B_value*twonorm*twonorm; + +// itk::DiffusionTensor3D< double > S; +// S[0] = g[0]*g[0]; +// S[1] = g[1]*g[0]; +// S[2] = g[2]*g[0]; +// S[3] = g[1]*g[1]; +// S[4] = g[2]*g[1]; +// S[5] = g[2]*g[2]; + +// double approx = 0; +// for (int i=0; iSetPixel(oit.GetIndex(), dPix); + + ++git; + } + + std::cout << "One Thread finished calculation" << std::endl; +} + +template< class TInPixelType, class TOutPixelType > +void +MultiTensorImageFilter< TInPixelType, TOutPixelType> +::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os,indent); +} + +} + +#endif // __itkMultiTensorImageFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkAdcImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkAdcImageFilter.h index 73a62bcaa5..87bbb5deff 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkAdcImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkAdcImageFilter.h @@ -1,135 +1,136 @@ /*=================================================================== 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 #include namespace itk{ /** \class AdcImageFilter */ template< class TInPixelType, class TOutPixelType > class AdcImageFilter : - public ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > + 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::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; - typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GradientContainerType; - typedef itk::Image< unsigned char, 3> ItkUcharImageType; - - /** 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( MaskImage, ItkUcharImageType::Pointer ) - itkSetMacro( FitSignal, bool ) - 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); - - bool m_FitSignal; - double m_B_value; - vnl_vector m_B_values; - vnl_vector m_Nonzero_B_values; - GradientContainerType m_GradientDirections; - ItkUcharImageType::Pointer m_MaskImage; - - double FitSingleVoxel( const typename InputImageType::PixelType &input); - - /** - * \brief The lestSquaresFunction struct for Non-Linear-Least-Squres fit of monoexponential model Si = S0*exp(-b*ADC) - */ - struct adcLeastSquaresFunction: public vnl_least_squares_function + typedef AdcImageFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< VectorImage< TInPixelType, 3 >, Image< TOutPixelType, 3 > > Superclass; + typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; + typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType::Pointer GradientContainerType; + typedef itk::Image< unsigned char, 3> ItkUcharImageType; + + /** 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( MaskImage, ItkUcharImageType::Pointer ) + itkSetMacro( FitSignal, bool ) + itkSetMacro( B_value, double ) + itkSetMacro( GradientDirections, GradientContainerType ) + + /** + * \brief The lestSquaresFunction struct for Non-Linear-Least-Squres fit of monoexponential model Si = S0*exp(-b*ADC) + */ + struct adcLeastSquaresFunction: public vnl_least_squares_function + { + + void set_S0(double val) { + S0 = val; + } - void set_S0(double val) - { - S0 = val; - } + void set_measurements(const vnl_vector& m) + { + measurements.set_size(m.size()); + measurements.copy_in(m.data_block()); + } - void set_measurements(const vnl_vector& m) - { - measurements.set_size(m.size()); - measurements.copy_in(m.data_block()); - } + void set_bvalues(const vnl_vector& x) + { + bValues.set_size(x.size()); + bValues.copy_in(x.data_block()); + } - void set_bvalues(const vnl_vector& x) - { - bValues.set_size(x.size()); - bValues.copy_in(x.data_block()); - } + vnl_vector measurements; + vnl_vector bValues; + double S0; + + adcLeastSquaresFunction(unsigned int number_of_measurements=1) : + vnl_least_squares_function(1, number_of_measurements, no_gradient) + { + } + + void f(const vnl_vector& x, vnl_vector& fx) override { - vnl_vector measurements; - vnl_vector bValues; - double S0; + const double & ADC = x[0]; - adcLeastSquaresFunction(unsigned int number_of_measurements=1) : - vnl_least_squares_function(1, number_of_measurements, no_gradient) + for(unsigned int s=0; s& x, vnl_vector& fx) override { + protected: + AdcImageFilter(); + ~AdcImageFilter() {} + void PrintSelf(std::ostream& os, Indent indent) const; - const double & ADC = x[0]; + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); + + bool m_FitSignal; + double m_B_value; + vnl_vector m_B_values; + vnl_vector m_Nonzero_B_values; + GradientContainerType m_GradientDirections; + ItkUcharImageType::Pointer m_MaskImage; + + double FitSingleVoxel( const typename InputImageType::PixelType &input); - for(unsigned int s=0; s #include #include #include #include #include #include #include #include namespace itk{ /** \brief Extracts principal eigenvectors of the input tensors */ template< class TTensorPixelType> class DiffusionTensorPrincipalDirectionImageFilter : public ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< unsigned char, 3 > > { public: typedef DiffusionTensorPrincipalDirectionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< unsigned char, 3 > > Superclass; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(DiffusionTensorPrincipalDirectionImageFilter, ImageToImageFilter) typedef TTensorPixelType TensorComponentType; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef itk::Image ItkUcharImgType; typedef vnl_vector_fixed< double, 3 > DirectionType; typedef Image< float, 4 > PeakImageType; void SetImage( const InputImageType *image ); // input itkSetMacro( MaskImage, ItkUcharImgType::Pointer) itkSetMacro( NormalizeVectors, bool) itkSetMacro( UsePolarCoordinates, bool) + itkSetMacro( FaThreshold, float) // output itkGetMacro( OutputFiberBundle, mitk::FiberBundle::Pointer) itkGetMacro( PeakImage, PeakImageType::Pointer) protected: DiffusionTensorPrincipalDirectionImageFilter(); ~DiffusionTensorPrincipalDirectionImageFilter() {} void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType ); void AfterThreadedGenerateData(); private: bool m_NormalizeVectors; ///< Normalizes the output vector to length 1 mitk::FiberBundle::Pointer m_OutputFiberBundle; ///< Vector field representation of the output vectors ItkUcharImgType::Pointer m_MaskImage; ///< Extraction is only performed inside of the binary mask PeakImageType::Pointer m_PeakImage; - float m_MaxEigenvalue; bool m_UsePolarCoordinates; + float m_FaThreshold; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionTensorPrincipalDirectionImageFilter.txx" #endif #endif //__itkDiffusionTensorPrincipalDirectionImageFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx index 3d1a4a771c..1fda8f5ff8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkDiffusionTensorPrincipalDirectionImageFilter.txx @@ -1,280 +1,282 @@ /*=================================================================== 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 __itkDiffusionTensorPrincipalDirectionImageFilter_txx #define __itkDiffusionTensorPrincipalDirectionImageFilter_txx #include #include #include #include "itkDiffusionTensorPrincipalDirectionImageFilter.h" #include "itkImageRegionConstIterator.h" #include "itkImageRegionConstIteratorWithIndex.h" #include "itkImageRegionIterator.h" #include "itkArray.h" #include "vnl/vnl_vector.h" #include #include #include #include #include #include #define _USE_MATH_DEFINES #include namespace itk { //#define QBALL_RECON_PI M_PI template< class TTensorPixelType> DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType>::DiffusionTensorPrincipalDirectionImageFilter() : m_NormalizeVectors(true) - , m_MaxEigenvalue(0.0) , m_UsePolarCoordinates(false) + , m_FaThreshold(0.2) { this->SetNumberOfRequiredInputs( 1 ); } template< class TTensorPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType>::BeforeThreadedGenerateData() { typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); Vector spacing3 = inputImagePointer->GetSpacing(); mitk::Point3D origin3 = inputImagePointer->GetOrigin(); itk::Matrix direction3 = inputImagePointer->GetDirection(); ImageRegion<3> imageRegion3 = inputImagePointer->GetLargestPossibleRegion(); if (m_MaskImage.IsNull()) { m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( spacing3 ); m_MaskImage->SetOrigin( origin3 ); m_MaskImage->SetDirection( direction3 ); m_MaskImage->SetRegions( imageRegion3 ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } typename OutputImageType::Pointer outputImage = OutputImageType::New(); outputImage->SetSpacing( spacing3 ); outputImage->SetOrigin( origin3 ); outputImage->SetDirection( direction3 ); outputImage->SetRegions( imageRegion3 ); outputImage->Allocate(); outputImage->FillBuffer(0); this->SetNthOutput(0, outputImage); itk::Vector spacing4; itk::Point origin4; itk::Matrix direction4; itk::ImageRegion<4> imageRegion4; spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin3[3] = 0; for (int r=0; r<3; r++) for (int c=0; c<3; c++) direction4[r][c] = direction3[r][c]; direction4[3][3] = 1; imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); imageRegion4.SetSize(3, 3); m_PeakImage = PeakImageType::New(); m_PeakImage->SetSpacing( spacing4 ); m_PeakImage->SetOrigin( origin4 ); m_PeakImage->SetDirection( direction4 ); m_PeakImage->SetRegions( imageRegion4 ); m_PeakImage->Allocate(); m_PeakImage->FillBuffer(0.0); } template< class TTensorPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType> ::AfterThreadedGenerateData() { vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); typename OutputImageType::Pointer numDirImage = static_cast< OutputImageType* >( this->ProcessObject::GetPrimaryOutput() ); ImageRegionConstIterator< OutputImageType > it(numDirImage, numDirImage->GetLargestPossibleRegion() ); mitk::Vector3D spacing = numDirImage->GetSpacing(); double minSpacing = spacing[0]; if (spacing[1]GetPixel(index)==0) { ++it; continue; } typename PeakImageType::IndexType peakIndex; peakIndex[0] = it.GetIndex()[0]; peakIndex[1] = it.GetIndex()[1]; peakIndex[2] = it.GetIndex()[2]; DirectionType dir; peakIndex[3] = 0; dir[0] = m_PeakImage->GetPixel(peakIndex); peakIndex[3] = 1; dir[1] = m_PeakImage->GetPixel(peakIndex); peakIndex[3] = 2; dir[2] = m_PeakImage->GetPixel(peakIndex); - if (!m_NormalizeVectors && m_MaxEigenvalue>0 && !m_UsePolarCoordinates) - dir /= m_MaxEigenvalue; - vtkSmartPointer container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = index[0]; center[1] = index[1]; center[2] = index[2]; itk::Point worldCenter; numDirImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); itk::Point worldStart; worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing; worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing; worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing; worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing; worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); ++it; } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData); } template< class TTensorPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType> ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType ) { typedef itk::DiffusionTensor3D TensorType; typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetPrimaryOutput()); ImageRegionIterator< OutputImageType > numDirectionsIterator(outputImage, outputRegionForThread); InputIteratorType tensorIterator(inputImagePointer, outputRegionForThread ); while( !tensorIterator.IsAtEnd() ) { typename InputImageType::IndexType index = tensorIterator.GetIndex(); if (m_MaskImage->GetPixel(index)==0) { ++tensorIterator; ++numDirectionsIterator; continue; } typename InputImageType::PixelType b = tensorIterator.Get(); TensorType tensor = b.GetDataPointer(); typename PeakImageType::IndexType peakIndex; peakIndex[0] = tensorIterator.GetIndex()[0]; peakIndex[1] = tensorIterator.GetIndex()[1]; peakIndex[2] = tensorIterator.GetIndex()[2]; typename TensorType::EigenValuesArrayType eigenvalues; typename TensorType::EigenVectorsMatrixType eigenvectors; if(tensor.GetTrace()!=0) { tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); vnl_vector_fixed vec; vec[0] = eigenvectors(2,0); vec[1] = eigenvectors(2,1); vec[2] = eigenvectors(2,2); + vec.normalize(); - if (!m_NormalizeVectors) - vec *= eigenvalues[2]; - - if (eigenvalues[2]>m_MaxEigenvalue) - m_MaxEigenvalue = eigenvalues[2]; + vnl_vector_fixed out; out.fill(0); - vnl_vector_fixed out; - if (m_UsePolarCoordinates) + float fa = tensor.GetFractionalAnisotropy(); + if (faM_PI) + if (m_UsePolarCoordinates) + { + if(vec[0] || vec[1] || vec[2]) { - out[1] = out[1] - M_PI; + out[0] = sqrt( vec[0] * vec[0] + vec[1] * vec[1] + vec[2] * vec[2] ); + out[1] = atan2( vec[1], vec[0] ); + out[2] = 0.5*M_PI - atan( vec[2] / sqrt( vec[0] * vec[0] + vec[1] * vec[1] ) ); + + if(out[1]>M_PI) + { + out[1] = out[1] - M_PI; + } + } + else + { + out[0] = 0; + out[1] = 0; + out[2] = 0; } } else { - out[0] = 0; - out[1] = 0; - out[2] = 0; + out = vec; } } - else - { - out = vec; - } peakIndex[3] = 0; m_PeakImage->SetPixel(peakIndex, out[0]); peakIndex[3] = 1; m_PeakImage->SetPixel(peakIndex, out[1]); peakIndex[3] = 2; m_PeakImage->SetPixel(peakIndex, out[2]); numDirectionsIterator.Set( 1 ); } ++numDirectionsIterator; ++tensorIterator; } std::cout << "One Thread finished extraction" << std::endl; } template< class TTensorPixelType> void DiffusionTensorPrincipalDirectionImageFilter< TTensorPixelType> ::PrintSelf(std::ostream& , Indent ) const { } } #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkDwiPhantomGenerationFilter.h b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkDwiPhantomGenerationFilter.h deleted file mode 100644 index b8b8b94771..0000000000 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkDwiPhantomGenerationFilter.h +++ /dev/null @@ -1,153 +0,0 @@ -/*=================================================================== - -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: Insight Segmentation & Registration Toolkit - Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $ - Language: C++ - Date: $Date: 2006-03-27 17:01:06 $ - Version: $Revision: 1.12 $ - - Copyright (c) Insight Software Consortium. All rights reserved. - See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm 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 __itkDwiPhantomGenerationFilter_h_ -#define __itkDwiPhantomGenerationFilter_h_ - -#include -#include -#include -#include -#include - -namespace itk{ - -/** -* \brief Generation of synthetic diffusion weighted images using a second rank tensor model. */ - -template< class TOutputScalarType > -class DwiPhantomGenerationFilter : public ImageSource< itk::VectorImage > -{ - -public: - - typedef DwiPhantomGenerationFilter Self; - typedef SmartPointer Pointer; - typedef SmartPointer ConstPointer; - typedef ImageSource< itk::VectorImage > Superclass; - - /** Method for creation through the object factory. */ - itkFactorylessNewMacro(Self) - itkCloneMacro(Self) - - /** Runtime information support. */ - itkTypeMacro(DwiPhantomGenerationFilter, ImageSource) - - typedef typename Superclass::OutputImageType OutputImageType; - typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - typedef itk::VectorContainer< int, double > AnglesContainerType; - typedef Vector GradientType; - typedef std::vector GradientListType; - typedef itk::Matrix MatrixType; - typedef itk::Image ItkUcharImgType; - typedef itk::Image ItkFloatImgType; - - typedef Image< Vector< float, 3 >, 3> ItkDirectionImage; - typedef VectorContainer< unsigned int, ItkDirectionImage::Pointer > ItkDirectionImageContainer; - - void SetGradientList(GradientListType gradientList){m_GradientList = gradientList;} - void SetSignalRegions(std::vector< ItkUcharImgType::Pointer > signalRegions){m_SignalRegions = signalRegions;} - void SetTensorFA(std::vector< float > faList){m_TensorFA = faList;} - void SetTensorADC(std::vector< float > adcList){m_TensorADC = adcList;} - void SetTensorWeight(std::vector< float > weightList){m_TensorWeight = weightList;} - void SetTensorDirection(std::vector< vnl_vector_fixed > directionList){m_TensorDirection = directionList;} - - // input parameters - itkSetMacro( BValue, float ) ///< signal parameter - itkSetMacro( SignalScale, float ) ///< scaling factor for signal - itkSetMacro( NoiseVariance, double ) ///< variance of rician noise - itkSetMacro( GreyMatterAdc, float ) ///< ADC of isotropic diffusion tensor - itkSetMacro( Spacing, mitk::Vector3D ) ///< parameter of output image - itkSetMacro( Origin, mitk::Point3D ) ///< parameter of output image - itkSetMacro( DirectionMatrix, MatrixType ) ///< parameter of output image - itkSetMacro( ImageRegion, ImageRegion<3> ) ///< parameter of output image - itkSetMacro( SimulateBaseline, bool ) ///< generate baseline image values as the l2 norm of the corresponding tensor used for the diffusion signal generation - - // output - itkGetMacro( DirectionImageContainer, ItkDirectionImageContainer::Pointer) ///< contains one vectorimage for each input ROI - itkGetMacro( NumDirectionsImage, ItkUcharImgType::Pointer) ///< contains number of directions per voxel - itkGetMacro( SNRImage, ItkFloatImgType::Pointer) ///< contains local SNR values - itkGetMacro( OutputFiberBundle, mitk::FiberBundle::Pointer) ///< output vector field - -protected: - DwiPhantomGenerationFilter(); - ~DwiPhantomGenerationFilter(){} - - void GenerateData(); - -private: - - // image variables - itk::Vector m_Spacing; - mitk::Point3D m_Origin; - itk::Matrix m_DirectionMatrix; - ImageRegion<3> m_ImageRegion; - ItkDirectionImageContainer::Pointer m_DirectionImageContainer; ///< contains one vectorimage for each input ROI - ItkUcharImgType::Pointer m_NumDirectionsImage; ///< contains number of directions per voxel - ItkFloatImgType::Pointer m_SNRImage; ///< contains local SNR values - mitk::FiberBundle::Pointer m_OutputFiberBundle; ///< output vector field - - // signal regions - std::vector< ItkUcharImgType::Pointer > m_SignalRegions; ///< binary images defining the regions for the signal generation - std::vector< float > m_TensorFA; ///< kernel tensor parameter - std::vector< float > m_TensorADC; ///< kernel tensor parameter - std::vector< float > m_TensorWeight; ///< weight factor of the signal - std::vector< vnl_vector_fixed > m_TensorDirection; ///< principal direction of the kernel tensor in the different ROIs - - // signal related variable - GradientListType m_GradientList; ///< list of used diffusion gradient directions - std::vector< itk::DiffusionTensor3D > m_TensorList; ///< kernel tensor rotated in the different directions - float m_BValue; ///< signal parameter - float m_SignalScale; ///< scaling factor for signal - Statistics::MersenneTwisterRandomVariateGenerator::Pointer m_RandGen; - int m_BaselineImages; ///< number of simulated baseline images - double m_MaxBaseline; ///< maximum value of the baseline image - double m_MeanBaseline; ///< mean value of the baseline image - double m_NoiseVariance; ///< variance of rician noise - bool m_AddNoiseFlag; //* Checkbox flag to switch noise on and off. */ - float m_GreyMatterAdc; ///< ADC of isotropic diffusion tensor - bool m_SimulateBaseline; ///< generate baseline image values as the l2 norm of the corresponding tensor used for the diffusion signal generation - TOutputScalarType m_DefaultBaseline; ///< default value for baseline image - - double GetTensorL2Norm(itk::DiffusionTensor3D& T); - void GenerateTensors(); ///< rotates kernel tensor in the direction of the input direction vectors - typename OutputImageType::PixelType SimulateMeasurement(itk::DiffusionTensor3D& tensor, float weight); - void AddNoise(typename OutputImageType::PixelType& pix); ///< Adds rician noise to the input pixel -}; - -} - -#ifndef ITK_MANUAL_INSTANTIATION -#include "itkDwiPhantomGenerationFilter.cpp" -#endif - -#endif //__itkDwiPhantomGenerationFilter_h_ - diff --git a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp index 1be4931f57..369f691042 100644 --- a/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/include/Algorithms/itkFiniteDiffOdfMaximaExtractionFilter.cpp @@ -1,548 +1,548 @@ /*=================================================================== 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 __itkFiniteDiffOdfMaximaExtractionFilter_cpp #define __itkFiniteDiffOdfMaximaExtractionFilter_cpp #include "itkFiniteDiffOdfMaximaExtractionFilter.h" #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace boost::math; namespace itk { static bool CompareVectors(const vnl_vector_fixed< double, 3 >& v1, const vnl_vector_fixed< double, 3 >& v2) { return (v1.magnitude()>v2.magnitude()); } template< class PixelType, int ShOrder, int NrOdfDirections > FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::FiniteDiffOdfMaximaExtractionFilter() : m_NormalizationMethod(MAX_VEC_NORM) , m_MaxNumPeaks(2) , m_PeakThreshold(0.4) , m_AbsolutePeakThreshold(0) , m_ClusteringThreshold(0.9) , m_AngularThreshold(0.7) , m_NumCoeffs((ShOrder*ShOrder + ShOrder + 2)/2 + ShOrder) , m_Toolkit(FSL) , m_ApplyDirectionMatrix(false) { this->SetNumberOfRequiredInputs(1); } template< class PixelType, int ShOrder, int NrOdfDirections > void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::FindCandidatePeaks(OdfType& odf, double thr, std::vector< DirectionType >& container) { double gfa = odf.GetGeneralizedFractionalAnisotropy(); //Find the peaks using a finite difference method bool flag = true; vnl_vector_fixed< bool, NrOdfDirections > used; used.fill(false); //Find the peaks for (int i=0; ithr && val*gfa>m_AbsolutePeakThreshold) // limit to one hemisphere ??? { flag = true; std::vector< int > neighbours = odf.GetNeighbors(i); for (unsigned int j=0; j std::vector< vnl_vector_fixed< double, 3 > > FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections>::MeanShiftClustering(std::vector< DirectionType >& inDirs) { std::vector< DirectionType > outDirs; if (inDirs.empty()) return inDirs; DirectionType oldMean, currentMean, workingMean; std::vector< int > touched; // initialize touched.resize(inDirs.size(), 0); bool free = true; currentMean = inDirs[0]; // initialize first seed while (free) { oldMean.fill(0.0); // start mean-shift clustering float angle = 0.0; int counter = 0; while ((currentMean-oldMean).magnitude()>0.0001) { counter = 0; oldMean = currentMean; workingMean = oldMean; workingMean.normalize(); currentMean.fill(0.0); for (unsigned int i=0; i=m_ClusteringThreshold) { currentMean += inDirs[i]; touched[i] = 1; counter++; } else if (-angle>=m_ClusteringThreshold) { currentMean -= inDirs[i]; touched[i] = 1; counter++; } } } // found stable mean if (counter>0) { float mag = currentMean.magnitude(); if (mag>0) { currentMean /= mag; outDirs.push_back(currentMean); } } // find next unused seed free = false; for (unsigned int i=0; i void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::BeforeThreadedGenerateData() { typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); itk::Vector spacing = ShCoeffImage->GetSpacing(); double minSpacing = spacing[0]; if (spacing[1]GetOrigin(); itk::Matrix direction = ShCoeffImage->GetDirection(); ImageRegion<3> imageRegion = ShCoeffImage->GetLargestPossibleRegion(); if (m_MaskImage.IsNotNull()) { origin = m_MaskImage->GetOrigin(); direction = m_MaskImage->GetDirection(); imageRegion = m_MaskImage->GetLargestPossibleRegion(); } itk::Vector spacing3 = ShCoeffImage->GetSpacing(); itk::Point origin3 = ShCoeffImage->GetOrigin(); itk::Matrix direction3 = ShCoeffImage->GetDirection(); itk::ImageRegion<3> imageRegion3 = ShCoeffImage->GetLargestPossibleRegion(); itk::Vector spacing4; itk::Point origin4; itk::Matrix direction4; itk::ImageRegion<4> imageRegion4; spacing4[0] = spacing3[0]; spacing4[1] = spacing3[1]; spacing4[2] = spacing3[2]; spacing4[3] = 1; - origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin[3] = 0; + origin4[0] = origin3[0]; origin4[1] = origin3[1]; origin4[2] = origin3[2]; origin4[3] = 0; for (int r=0; r<3; r++) for (int c=0; c<3; c++) direction4[r][c] = direction3[r][c]; direction4[3][3] = 1; imageRegion4.SetSize(0, imageRegion3.GetSize()[0]); imageRegion4.SetSize(1, imageRegion3.GetSize()[1]); imageRegion4.SetSize(2, imageRegion3.GetSize()[2]); imageRegion4.SetSize(3, m_MaxNumPeaks*3); m_PeakImage = PeakImageType::New(); m_PeakImage->SetSpacing( spacing4 ); m_PeakImage->SetOrigin( origin4 ); m_PeakImage->SetDirection( direction4 ); m_PeakImage->SetRegions( imageRegion4 ); m_PeakImage->Allocate(); m_PeakImage->FillBuffer(0.0); if (m_MaskImage.IsNull()) { m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( spacing ); m_MaskImage->SetOrigin( origin ); m_MaskImage->SetDirection( direction ); m_MaskImage->SetRegions( imageRegion ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_NumDirectionsImage = ItkUcharImgType::New(); m_NumDirectionsImage->SetSpacing( spacing ); m_NumDirectionsImage->SetOrigin( origin ); m_NumDirectionsImage->SetDirection( direction ); m_NumDirectionsImage->SetRegions( imageRegion ); m_NumDirectionsImage->Allocate(); m_NumDirectionsImage->FillBuffer(0); // calculate SH basis OdfType odf; vnl_matrix< double > sphCoords; std::vector< DirectionType > dirs; for (int i=0; i void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::AfterThreadedGenerateData() { MITK_INFO << "Generating vector field"; vtkSmartPointer m_VtkCellArray = vtkSmartPointer::New(); vtkSmartPointer m_VtkPoints = vtkSmartPointer::New(); typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, ShCoeffImage->GetLargestPossibleRegion() ); mitk::Vector3D spacing = ShCoeffImage->GetSpacing(); double minSpacing = spacing[0]; if (spacing[1]GetLargestPossibleRegion().GetSize()[0]*ShCoeffImage->GetLargestPossibleRegion().GetSize()[1]*ShCoeffImage->GetLargestPossibleRegion().GetSize()[2]; boost::progress_display disp(maxProgress); while( !cit.IsAtEnd() ) { ++disp; typename CoefficientImageType::IndexType idx3 = cit.GetIndex(); if (m_MaskImage->GetPixel(idx3)==0) { ++cit; continue; } itk::Index<4> idx4; idx4[0] = idx3[0]; idx4[1] = idx3[1]; idx4[2] = idx3[2]; for (unsigned int i=0; iGetPixel(idx4); idx4[3] = i*3 + 1; dir[1] = m_PeakImage->GetPixel(idx4); idx4[3] = i*3 + 2; dir[2] = m_PeakImage->GetPixel(idx4); vtkSmartPointer container = vtkSmartPointer::New(); itk::ContinuousIndex center; center[0] = idx3[0]; center[1] = idx3[1]; center[2] = idx3[2]; itk::Point worldCenter; m_MaskImage->TransformContinuousIndexToPhysicalPoint( center, worldCenter ); itk::Point worldStart; worldStart[0] = worldCenter[0]-dir[0]/2 * minSpacing; worldStart[1] = worldCenter[1]-dir[1]/2 * minSpacing; worldStart[2] = worldCenter[2]-dir[2]/2 * minSpacing; vtkIdType id = m_VtkPoints->InsertNextPoint(worldStart.GetDataPointer()); container->GetPointIds()->InsertNextId(id); itk::Point worldEnd; worldEnd[0] = worldCenter[0]+dir[0]/2 * minSpacing; worldEnd[1] = worldCenter[1]+dir[1]/2 * minSpacing; worldEnd[2] = worldCenter[2]+dir[2]/2 * minSpacing; id = m_VtkPoints->InsertNextPoint(worldEnd.GetDataPointer()); container->GetPointIds()->InsertNextId(id); m_VtkCellArray->InsertNextCell(container); } ++cit; } vtkSmartPointer directionsPolyData = vtkSmartPointer::New(); directionsPolyData->SetPoints(m_VtkPoints); directionsPolyData->SetLines(m_VtkCellArray); m_OutputFiberBundle = mitk::FiberBundle::New(directionsPolyData); } template< class PixelType, int ShOrder, int NrOdfDirections > void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::ThreadedGenerateData( const OutputImageRegionType& outputRegionForThread, ThreadIdType threadID ) { typename CoefficientImageType::Pointer ShCoeffImage = static_cast< CoefficientImageType* >( this->ProcessObject::GetInput(0) ); ImageRegionConstIterator< CoefficientImageType > cit(ShCoeffImage, outputRegionForThread ); OdfType odf; while( !cit.IsAtEnd() ) { typename CoefficientImageType::IndexType idx3 = cit.GetIndex(); if (m_MaskImage->GetPixel(idx3)==0) { ++cit; continue; } CoefficientPixelType c = cit.Get(); // calculate ODF double max = 0; odf.Fill(0.0); for (int i=0; imax) max = odf[i]; } if (max<0.0001) { ++cit; continue; } std::vector< DirectionType > candidates, peaks, temp; peaks.clear(); max *= m_PeakThreshold; // relative threshold FindCandidatePeaks(odf, max, candidates); // find all local maxima candidates = MeanShiftClustering(candidates); // cluster maxima vnl_matrix< double > shBasis, sphCoords; Cart2Sph(candidates, sphCoords); // convert candidate peaks to spherical angles shBasis = CalcShBasis(sphCoords); // evaluate spherical harmonics at each peak max = 0.0; for (unsigned int i=0; imax) max = val; peaks.push_back(candidates[i]*val); } std::sort( peaks.begin(), peaks.end(), CompareVectors ); // sort peaks // kick out directions to close to a larger direction (too far away to cluster but too close to keep) unsigned int m = peaks.size(); if ( m>m_MaxNumPeaks ) m = m_MaxNumPeaks; for (unsigned int i=0; im_AngularThreshold && val idx4; idx4[0] = idx3[0]; idx4[1] = idx3[1]; idx4[2] = idx3[2]; // fill output image unsigned int num = peaks.size(); if ( num>m_MaxNumPeaks ) num = m_MaxNumPeaks; for (unsigned int i=0; iGetDirection()*dir; if (m_FlipX) dir[0] = -dir[0]; if (m_FlipY) dir[1] = -dir[1]; if (m_FlipZ) dir[2] = -dir[2]; for (unsigned int j = 0; j<3; j++) { idx4[3] = i*3 + j; m_PeakImage->SetPixel(idx4, dir[j]); } } m_NumDirectionsImage->SetPixel(idx3, num); ++cit; } MITK_INFO << "Thread " << threadID << " finished extraction"; } // convert cartesian to spherical coordinates template< class PixelType, int ShOrder, int NrOdfDirections > void FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::Cart2Sph(const std::vector< DirectionType >& dir, vnl_matrix& sphCoords) { sphCoords.set_size(dir.size(), 2); for (unsigned int i=0; i vnl_matrix FiniteDiffOdfMaximaExtractionFilter< PixelType, ShOrder, NrOdfDirections> ::CalcShBasis(vnl_matrix& sphCoords) { int M = sphCoords.rows(); int j, m; double mag, plm; vnl_matrix shBasis; shBasis.set_size(M, m_NumCoeffs); for (int p=0; p(l,abs(m),cos(sphCoords(p,0))); mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; if (m<0) shBasis(p,j) = sqrt(2.0)*mag*cos(fabs((double)m)*sphCoords(p,1)); else if (m==0) shBasis(p,j) = mag; else shBasis(p,j) = pow(-1.0, m)*sqrt(2.0)*mag*sin(m*sphCoords(p,1)); break; case MRTRIX: plm = legendre_p(l,abs(m),-cos(sphCoords(p,0))); mag = sqrt((double)(2*l+1)/(4.0*M_PI)*factorial(l-abs(m))/factorial(l+abs(m)))*plm; if (m>0) shBasis(p,j) = mag*cos(m*sphCoords(p,1)); else if (m==0) shBasis(p,j) = mag; else shBasis(p,j) = mag*sin(-m*sphCoords(p,1)); break; } j++; } } return shBasis; } } #endif // __itkFiniteDiffOdfMaximaExtractionFilter_cpp diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp index 94cf74abc2..2c5350d36a 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.odfpeaks/src/internal/QmitkOdfMaximaExtractionView.cpp @@ -1,477 +1,478 @@ /*=================================================================== 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. ===================================================================*/ //misc #define _USE_MATH_DEFINES #include #include // Blueberry #include #include // Qmitk #include "QmitkOdfMaximaExtractionView.h" // MITK #include #include #include #include #include // ITK #include #include #include #include #include #include #include #include #include #include #include #include // Qt #include const std::string QmitkOdfMaximaExtractionView::VIEW_ID = "org.mitk.views.odfmaximaextractionview"; using namespace mitk; QmitkOdfMaximaExtractionView::QmitkOdfMaximaExtractionView() : m_Controls(nullptr) { } // Destructor QmitkOdfMaximaExtractionView::~QmitkOdfMaximaExtractionView() { } void QmitkOdfMaximaExtractionView::CreateQtPartControl(QWidget *parent) { // build up qt view, unless already done if (!m_Controls) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkOdfMaximaExtractionViewControls; m_Controls->setupUi(parent); connect((QObject*)m_Controls->m_StartPeakExtractionButton, SIGNAL(clicked()), (QObject*) this, SLOT(StartPeakExtraction())); connect((QObject*)m_Controls->m_ImportShCoeffs, SIGNAL(clicked()), (QObject*) this, SLOT(ConvertShCoeffs())); m_Controls->m_MaskBox->SetDataStorage(this->GetDataStorage()); m_Controls->m_ImageBox->SetDataStorage(this->GetDataStorage()); mitk::TNodePredicateDataType::Pointer isMitkImage = mitk::TNodePredicateDataType::New(); mitk::NodePredicateNot::Pointer isDwi = mitk::NodePredicateNot::New(mitk::NodePredicateIsDWI::New()); mitk::NodePredicateNot::Pointer isQbi = mitk::NodePredicateNot::New(mitk::NodePredicateDataType::New("QBallImage")); mitk::NodePredicateAnd::Pointer unwanted = mitk::NodePredicateAnd::New(isQbi, isDwi); mitk::NodePredicateDimension::Pointer dim3 = mitk::NodePredicateDimension::New(3); mitk::NodePredicateProperty::Pointer isBinaryPredicate = mitk::NodePredicateProperty::New("binary", mitk::BoolProperty::New(true)); m_Controls->m_MaskBox->SetPredicate(mitk::NodePredicateAnd::New(mitk::NodePredicateAnd::New(unwanted, dim3), isBinaryPredicate)); m_Controls->m_ImageBox->SetPredicate(mitk::NodePredicateAnd::New(mitk::NodePredicateAnd::New(unwanted, isMitkImage), mitk::NodePredicateNot::New(isBinaryPredicate))); m_Controls->m_MaskBox->SetZeroEntryText("--"); m_Controls->m_ImageBox->SetZeroEntryText("--"); connect( (QObject*)(m_Controls->m_ImageBox), SIGNAL(OnSelectionChanged(const mitk::DataNode*)), this, SLOT(OnImageSelectionChanged()) ); m_Controls->m_StartPeakExtractionButton->setVisible(false); m_Controls->m_ImportShCoeffs->setVisible(false); } } void QmitkOdfMaximaExtractionView::SetFocus() { } void QmitkOdfMaximaExtractionView::StartPeakExtraction() { if (dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()) != nullptr) { StartTensorPeakExtraction(dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData())); } else { StartMaximaExtraction(dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData())); } } template void QmitkOdfMaximaExtractionView::TemplatedConvertShCoeffs(mitk::Image* mitkImg) { typedef itk::ShCoefficientImageImporter< float, shOrder > FilterType; typedef mitk::ImageToItk< itk::Image< float, 4 > > CasterType; CasterType::Pointer caster = CasterType::New(); caster->SetInput(mitkImg); caster->Update(); typename FilterType::Pointer filter = FilterType::New(); switch (m_Controls->m_ToolkitBox->currentIndex()) { case 0: filter->SetToolkit(FilterType::FSL); break; case 1: filter->SetToolkit(FilterType::MRTRIX); break; default: filter->SetToolkit(FilterType::FSL); } filter->SetInputImage(caster->GetOutput()); filter->GenerateData(); typename FilterType::QballImageType::Pointer itkQbi = filter->GetQballImage(); typename FilterType::CoefficientImageType::Pointer itkCi = filter->GetCoefficientImage(); { mitk::Image::Pointer img = mitk::Image::New(); img->InitializeByItk(itkCi.GetPointer()); img->SetVolume(itkCi->GetBufferPointer()); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_ShCoefficientImage_Imported"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); } { mitk::QBallImage::Pointer img = mitk::QBallImage::New(); img->InitializeByItk(itkQbi.GetPointer()); img->SetVolume(itkQbi->GetBufferPointer()); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_OdfImage_Imported"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); } } void QmitkOdfMaximaExtractionView::ConvertShCoeffs() { if (m_Controls->m_ImageBox->GetSelectedNode().IsNull()) return; Image::Pointer mitkImg = dynamic_cast(m_Controls->m_ImageBox->GetSelectedNode()->GetData()); if (mitkImg->GetDimension() != 4 && mitkImg->GetLargestPossibleRegion().GetSize()[3]<6) { MITK_INFO << "wrong image type (need 4 dimensions)"; return; } int nrCoeffs = mitkImg->GetLargestPossibleRegion().GetSize()[3]; // // solve bx² + cx + d = 0 = shOrder² + 2*shOrder + 2-2*neededCoeffs; // int c = 3, d = 2 - 2 * nrCoeffs; // double D = c*c - 4 * d; // int shOrder; // if (D>0) // { // shOrder = (-c + sqrt(D)) / 2.0; // if (shOrder<0) // shOrder = (-c - sqrt(D)) / 2.0; // } // else if (D == 0) // shOrder = -c / 2.0; switch (nrCoeffs) { case 6: TemplatedConvertShCoeffs<2>(mitkImg); break; case 15: TemplatedConvertShCoeffs<4>(mitkImg); break; case 28: TemplatedConvertShCoeffs<6>(mitkImg); break; case 45: TemplatedConvertShCoeffs<8>(mitkImg); break; case 66: TemplatedConvertShCoeffs<10>(mitkImg); break; case 91: TemplatedConvertShCoeffs<12>(mitkImg); break; default : QMessageBox::warning(nullptr, "Error", "Only spherical harmonics orders 2-12 are supported.", QMessageBox::Ok); } } void QmitkOdfMaximaExtractionView::StartTensorPeakExtraction(mitk::TensorImage* img) { typedef itk::DiffusionTensorPrincipalDirectionImageFilter< float > MaximaExtractionFilterType; MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); filter->SetUsePolarCoordinates(false); mitk::BaseGeometry::Pointer geometry; try{ ItkTensorImage::Pointer itkImage = ItkTensorImage::New(); CastToItkImage(img, itkImage); filter->SetInput(itkImage); geometry = img->GetGeometry(); } catch (itk::ExceptionObject &e) { MITK_INFO << "wrong image type: " << e.what(); QMessageBox::warning(nullptr, "Wrong pixel type", "Could not perform Tensor Principal Direction Extraction due to Image has wrong pixel type.", QMessageBox::Ok); return; } if (m_Controls->m_MaskBox->GetSelectedNode().IsNotNull()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); Image::Pointer mitkMaskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); CastToItkImage(mitkMaskImg, itkMaskImage); filter->SetMaskImage(itkMaskImage); } if (m_Controls->m_NormalizationBox->currentIndex() == 0) filter->SetNormalizeVectors(false); + filter->SetFaThreshold(m_Controls->m_AbsoluteThresholdBox->value()); filter->Update(); if (m_Controls->m_OutputDirectionImagesBox->isChecked()) { MaximaExtractionFilterType::PeakImageType::Pointer itkImg = filter->GetPeakImage(); mitk::Image::Pointer img = mitk::Image::New(); CastToMitkImage(itkImg, img); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_PrincipalDirection"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); } if (m_Controls->m_OutputNumDirectionsBox->isChecked()) { ItkUcharImgType::Pointer numDirImage = filter->GetOutput(); mitk::Image::Pointer image2 = mitk::Image::New(); image2->InitializeByItk(numDirImage.GetPointer()); image2->SetVolume(numDirImage->GetBufferPointer()); DataNode::Pointer node2 = DataNode::New(); node2->SetData(image2); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_NumDirections"; node2->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node2, m_Controls->m_ImageBox->GetSelectedNode()); } if (m_Controls->m_OutputVectorFieldBox->isChecked()) { mitk::Vector3D outImageSpacing = geometry->GetSpacing(); float minSpacing = 1; if (outImageSpacing[0]GetOutputFiberBundle(); // directions->SetGeometry(geometry); DataNode::Pointer node = DataNode::New(); node->SetData(directions); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_VectorField"; node->SetName(name.toStdString().c_str()); node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(minSpacing)); node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); } } template void QmitkOdfMaximaExtractionView::StartMaximaExtraction(Image *image) { typedef itk::FiniteDiffOdfMaximaExtractionFilter< float, shOrder, 20242 > MaximaExtractionFilterType; typename MaximaExtractionFilterType::Pointer filter = MaximaExtractionFilterType::New(); switch (m_Controls->m_ToolkitBox->currentIndex()) { case 0: filter->SetToolkit(MaximaExtractionFilterType::FSL); break; case 1: filter->SetToolkit(MaximaExtractionFilterType::MRTRIX); break; default: filter->SetToolkit(MaximaExtractionFilterType::FSL); } mitk::BaseGeometry::Pointer geometry; try{ typedef ImageToItk< typename MaximaExtractionFilterType::CoefficientImageType > CasterType; typename CasterType::Pointer caster = CasterType::New(); caster->SetInput(image); caster->Update(); filter->SetInput(caster->GetOutput()); geometry = image->GetGeometry(); } catch (itk::ExceptionObject &e) { MITK_INFO << "wrong image type: " << e.what(); QMessageBox::warning(nullptr, "Wrong pixel type", "Could not perform Finite Differences Extraction due to Image has wrong pixel type.", QMessageBox::Ok); return; } filter->SetAngularThreshold(cos((float)m_Controls->m_AngularThreshold->value()*M_PI / 180)); filter->SetClusteringThreshold(cos((float)m_Controls->m_ClusteringAngleBox->value()*M_PI / 180)); filter->SetMaxNumPeaks(m_Controls->m_MaxNumPeaksBox->value()); filter->SetPeakThreshold(m_Controls->m_PeakThresholdBox->value()); filter->SetAbsolutePeakThreshold(m_Controls->m_AbsoluteThresholdBox->value()); if (m_Controls->m_MaskBox->GetSelectedNode().IsNotNull()) { ItkUcharImgType::Pointer itkMaskImage = ItkUcharImgType::New(); Image::Pointer mitkMaskImg = dynamic_cast(m_Controls->m_MaskBox->GetSelectedNode()->GetData()); CastToItkImage(mitkMaskImg, itkMaskImage); filter->SetMaskImage(itkMaskImage); } switch (m_Controls->m_NormalizationBox->currentIndex()) { case 0: filter->SetNormalizationMethod(MaximaExtractionFilterType::NO_NORM); break; case 1: filter->SetNormalizationMethod(MaximaExtractionFilterType::MAX_VEC_NORM); break; case 2: filter->SetNormalizationMethod(MaximaExtractionFilterType::SINGLE_VEC_NORM); break; } filter->Update(); if (m_Controls->m_OutputDirectionImagesBox->isChecked()) { typename MaximaExtractionFilterType::PeakImageType::Pointer itkImg = filter->GetPeakImage(); mitk::Image::Pointer img = mitk::Image::New(); CastToMitkImage(itkImg, img); DataNode::Pointer node = DataNode::New(); node->SetData(img); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_PEAKS"; node->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); } if (m_Controls->m_OutputNumDirectionsBox->isChecked()) { ItkUcharImgType::Pointer numDirImage = filter->GetNumDirectionsImage(); mitk::Image::Pointer image2 = mitk::Image::New(); CastToMitkImage(numDirImage, image2); DataNode::Pointer node2 = DataNode::New(); node2->SetData(image2); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_NUM_DIRECTIONS"; node2->SetName(name.toStdString().c_str()); GetDataStorage()->Add(node2, m_Controls->m_ImageBox->GetSelectedNode()); } if (m_Controls->m_OutputVectorFieldBox->isChecked()) { mitk::Vector3D outImageSpacing = geometry->GetSpacing(); float minSpacing = 1; if (outImageSpacing[0]GetOutputFiberBundle(); // directions->SetGeometry(geometry); DataNode::Pointer node = DataNode::New(); node->SetData(directions); QString name(m_Controls->m_ImageBox->GetSelectedNode()->GetName().c_str()); name += "_VECTOR_FIELD"; node->SetName(name.toStdString().c_str()); node->SetProperty("Fiber2DSliceThickness", mitk::FloatProperty::New(minSpacing)); node->SetProperty("Fiber2DfadeEFX", mitk::BoolProperty::New(false)); GetDataStorage()->Add(node, m_Controls->m_ImageBox->GetSelectedNode()); } } void QmitkOdfMaximaExtractionView::StartMaximaExtraction(Image* img) { mitk::PixelType pixT = img->GetPixelType(); switch (pixT.GetNumberOfComponents()) { case 6: StartMaximaExtraction<2>(img); break; case 15: StartMaximaExtraction<4>(img); break; case 28: StartMaximaExtraction<6>(img); break; case 45: StartMaximaExtraction<8>(img); break; case 66: StartMaximaExtraction<10>(img); break; case 91: StartMaximaExtraction<12>(img); break; default : QMessageBox::warning(nullptr, "Error", "Only spherical harmonics orders 2-12 are supported.", QMessageBox::Ok); } } void QmitkOdfMaximaExtractionView::OnSelectionChanged(berry::IWorkbenchPart::Pointer , const QList& nodes) { (void) nodes; this->OnImageSelectionChanged(); } void QmitkOdfMaximaExtractionView::OnImageSelectionChanged() { m_Controls->m_StartPeakExtractionButton->setVisible(false); m_Controls->m_ImportShCoeffs->setVisible(false); mitk::DataNode::Pointer node = m_Controls->m_ImageBox->GetSelectedNode(); if (node.IsNull()) return; Image::Pointer img = dynamic_cast(node->GetData()); if (img->GetDimension()==4) m_Controls->m_ImportShCoeffs->setVisible(true); else m_Controls->m_StartPeakExtractionButton->setVisible(true); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.cpp index e6d167c95d..9993bb4e8a 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.cpp @@ -1,632 +1,760 @@ /*=================================================================== 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 "QmitkDiffusionQuantificationView.h" #include "mitkDiffusionImagingConfigure.h" #include "itkTimeProbe.h" #include "itkImage.h" #include "mitkNodePredicateDataType.h" #include "mitkDataNodeObject.h" #include "mitkQBallImage.h" #include #include "mitkImageCast.h" #include "mitkStatusBar.h" #include "itkDiffusionQballGeneralizedFaImageFilter.h" #include "itkShiftScaleImageFilter.h" #include "itkTensorFractionalAnisotropyImageFilter.h" #include "itkTensorRelativeAnisotropyImageFilter.h" #include "itkTensorDerivedMeasurementsFilter.h" #include "QmitkDataStorageComboBox.h" #include #include "berryIWorkbenchWindow.h" #include "berryISelectionService.h" #include #include +#include +#include #include #include +#include +#include const std::string QmitkDiffusionQuantificationView::VIEW_ID = "org.mitk.views.diffusionquantification"; QmitkDiffusionQuantificationView::QmitkDiffusionQuantificationView() : QmitkAbstractView(), m_Controls(nullptr) { m_QBallImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); m_DwImages = mitk::DataStorage::SetOfObjects::New(); } QmitkDiffusionQuantificationView::~QmitkDiffusionQuantificationView() { } void QmitkDiffusionQuantificationView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkDiffusionQuantificationViewControls; m_Controls->setupUi(parent); this->CreateConnections(); GFACheckboxClicked(); #ifndef DIFFUSION_IMAGING_EXTENDED m_Controls->m_StandardGFACheckbox->setVisible(false); m_Controls->frame_3->setVisible(false); m_Controls->m_CurvatureButton->setVisible(false); #endif + + m_Controls->m_BallStickButton->setVisible(false); + m_Controls->m_MultiTensorButton->setVisible(false); } } void QmitkDiffusionQuantificationView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_StandardGFACheckbox), SIGNAL(clicked()), this, SLOT(GFACheckboxClicked()) ); connect( (QObject*)(m_Controls->m_GFAButton), SIGNAL(clicked()), this, SLOT(GFA()) ); connect( (QObject*)(m_Controls->m_CurvatureButton), SIGNAL(clicked()), this, SLOT(Curvature()) ); connect( (QObject*)(m_Controls->m_FAButton), SIGNAL(clicked()), this, SLOT(FA()) ); connect( (QObject*)(m_Controls->m_RAButton), SIGNAL(clicked()), this, SLOT(RA()) ); connect( (QObject*)(m_Controls->m_ADButton), SIGNAL(clicked()), this, SLOT(AD()) ); connect( (QObject*)(m_Controls->m_RDButton), SIGNAL(clicked()), this, SLOT(RD()) ); connect( (QObject*)(m_Controls->m_MDButton), SIGNAL(clicked()), this, SLOT(MD()) ); connect( (QObject*)(m_Controls->m_MdDwiButton), SIGNAL(clicked()), this, SLOT(MD_DWI()) ); connect( (QObject*)(m_Controls->m_AdcDwiButton), SIGNAL(clicked()), this, SLOT(ADC_DWI()) ); connect( (QObject*)(m_Controls->m_ClusteringAnisotropy), SIGNAL(clicked()), this, SLOT(ClusterAnisotropy()) ); + connect( (QObject*)(m_Controls->m_BallStickButton), SIGNAL(clicked()), this, SLOT(DoBallStickCalculation()) ); + connect( (QObject*)(m_Controls->m_MultiTensorButton), SIGNAL(clicked()), this, SLOT(DoMultiTensorCalculation()) ); } } void QmitkDiffusionQuantificationView::SetFocus() { m_Controls->m_ScaleImageValuesBox->setFocus(); } void QmitkDiffusionQuantificationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { m_QBallImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); m_DwImages = mitk::DataStorage::SetOfObjects::New(); bool foundQBIVolume = false; bool foundTensorVolume = false; bool foundDwVolume = false; mitk::DataNode::Pointer selNode = nullptr; int c=0; for (mitk::DataNode::Pointer node: nodes) { if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { foundQBIVolume = true; m_QBallImages->push_back(node); selNode = node; c++; } else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { foundTensorVolume = true; m_TensorImages->push_back(node); selNode = node; c++; } else if( node.IsNotNull() && dynamic_cast(node->GetData()) && mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(dynamic_cast(node->GetData())) ) { foundDwVolume = true; m_DwImages->push_back(node); selNode = node; c++; } } m_Controls->m_GFAButton->setEnabled(foundQBIVolume); m_Controls->m_CurvatureButton->setEnabled(foundQBIVolume); if (c>0) { m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_InputImageLabel->setText(selNode->GetName().c_str()); } else { m_Controls->m_InputData->setTitle("Please Select Input Data"); m_Controls->m_InputImageLabel->setText("mandatory"); } m_Controls->m_FAButton->setEnabled(foundTensorVolume); m_Controls->m_RAButton->setEnabled(foundTensorVolume); m_Controls->m_ADButton->setEnabled(foundTensorVolume); m_Controls->m_RDButton->setEnabled(foundTensorVolume); m_Controls->m_MDButton->setEnabled(foundTensorVolume); m_Controls->m_ClusteringAnisotropy->setEnabled(foundTensorVolume); m_Controls->m_AdcDwiButton->setEnabled(foundDwVolume); m_Controls->m_MdDwiButton->setEnabled(foundDwVolume); + m_Controls->m_BallStickButton->setEnabled(foundDwVolume); + m_Controls->m_MultiTensorButton->setEnabled(foundDwVolume); } void QmitkDiffusionQuantificationView::ADC_DWI() { DoAdcCalculation(true); } void QmitkDiffusionQuantificationView::MD_DWI() { DoAdcCalculation(false); } +void QmitkDiffusionQuantificationView::DoBallStickCalculation() +{ + mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_DwImages->end() ); + mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_DwImages->begin() ); + + while ( itemiter != itemiterend ) // for all items + { + mitk::DataNode* node = *itemiter; + mitk::Image::Pointer image = dynamic_cast(node->GetData()); + + typedef itk::BallAndSticksImageFilter< short, double > FilterType; + + ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); + mitk::CastToItkImage(image, itkVectorImagePointer); + FilterType::Pointer filter = FilterType::New(); + filter->SetInput( itkVectorImagePointer ); + + filter->SetGradientDirections( static_cast + ( image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ) + ->GetGradientDirectionsContainer() ); + + filter->SetB_value( static_cast + (image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() ) + ->GetValue() ); + + filter->Update(); + + mitk::Image::Pointer newImage = mitk::Image::New(); + newImage->InitializeByItk( filter->GetOutput() ); + newImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); + mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); + imageNode->SetData( newImage ); + QString name = node->GetName().c_str(); + + imageNode->SetName((name+"_f").toStdString().c_str()); + GetDataStorage()->Add(imageNode, node); + + { + FilterType::PeakImageType::Pointer itkImg = filter->GetPeakImage(); + mitk::Image::Pointer newImage = mitk::Image::New(); + CastToMitkImage(itkImg, newImage); + + mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); + imageNode->SetData( newImage ); + QString name = node->GetName().c_str(); + + imageNode->SetName((name+"_Sticks").toStdString().c_str()); + GetDataStorage()->Add(imageNode, node); + } + + { + mitk::Image::Pointer dOut = mitk::GrabItkImageMemory( filter->GetOutDwi().GetPointer() ); + + dOut->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ); + dOut->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() ); + mitk::DiffusionPropertyHelper propertyHelper( dOut ); + propertyHelper.InitializeImage(); + + mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); + imageNode->SetData( dOut ); + QString name = node->GetName().c_str(); + imageNode->SetName((name+"_Estimated-DWI").toStdString().c_str()); + GetDataStorage()->Add(imageNode, node); + } + + ++itemiter; + } +} + +void QmitkDiffusionQuantificationView::DoMultiTensorCalculation() +{ + mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_DwImages->end() ); + mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_DwImages->begin() ); + + while ( itemiter != itemiterend ) // for all items + { + mitk::DataNode* node = *itemiter; + mitk::Image::Pointer image = dynamic_cast(node->GetData()); + + typedef itk::MultiTensorImageFilter< short, double > FilterType; + + ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); + mitk::CastToItkImage(image, itkVectorImagePointer); + FilterType::Pointer filter = FilterType::New(); + filter->SetInput( itkVectorImagePointer ); + + filter->SetGradientDirections( static_cast + ( image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ) + ->GetGradientDirectionsContainer() ); + + filter->SetB_value( static_cast + (image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() ) + ->GetValue() ); + + filter->Update(); + + typedef itk::Image, 3> TensorImageType; + for (int i=0; iGetTensorImages().at(i); + mitk::TensorImage::Pointer image = mitk::TensorImage::New(); + image->InitializeByItk( tensorImage.GetPointer() ); + image->SetVolume( tensorImage->GetBufferPointer() ); + + mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); + imageNode->SetData( image ); + QString name = node->GetName().c_str(); + name.append("_Tensor"); + name.append(boost::lexical_cast(i).c_str()); + imageNode->SetName(name.toStdString().c_str()); + GetDataStorage()->Add(imageNode, node); + } + + ++itemiter; + } +} + void QmitkDiffusionQuantificationView::DoAdcCalculation(bool fit) { mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_DwImages->end() ); mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_DwImages->begin() ); while ( itemiter != itemiterend ) // for all items { mitk::DataNode* node = *itemiter; mitk::Image::Pointer image = dynamic_cast(node->GetData()); typedef itk::AdcImageFilter< short, double > FilterType; ItkDwiType::Pointer itkVectorImagePointer = ItkDwiType::New(); mitk::CastToItkImage(image, itkVectorImagePointer); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkVectorImagePointer ); filter->SetGradientDirections( static_cast ( image->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() ) ->GetGradientDirectionsContainer() ); filter->SetB_value( static_cast (image->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() ) ->GetValue() ); filter->SetFitSignal(fit); filter->Update(); typedef itk::ShiftScaleImageFilter::OutputImageType, itk::AdcImageFilter< short, double >::OutputImageType> ShiftScaleFilterType; ShiftScaleFilterType::Pointer multi = ShiftScaleFilterType::New(); multi->SetShift(0.0); multi->SetScale(m_Controls->m_ScaleImageValuesBox->value()); multi->SetInput(filter->GetOutput()); multi->Update(); mitk::Image::Pointer newImage = mitk::Image::New(); newImage->InitializeByItk( multi->GetOutput() ); newImage->SetVolume( multi->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( newImage ); QString name = node->GetName().c_str(); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType( mitk::LookupTable::JET_TRANSPARENT ); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable( lut ); imageNode->SetProperty("LookupTable", lut_prop ); if (fit) imageNode->SetName((name+"_ADC").toStdString().c_str()); else imageNode->SetName((name+"_MD").toStdString().c_str()); GetDataStorage()->Add(imageNode, node); ++itemiter; } } void QmitkDiffusionQuantificationView::GFACheckboxClicked() { m_Controls->frame_2->setVisible(m_Controls->m_StandardGFACheckbox->isChecked()); } void QmitkDiffusionQuantificationView::GFA() { if(m_Controls->m_StandardGFACheckbox->isChecked()) { QBIQuantify(13); } else { QBIQuantify(0); } } void QmitkDiffusionQuantificationView::Curvature() { QBIQuantify(12); } void QmitkDiffusionQuantificationView::FA() { TensorQuantify(0); } void QmitkDiffusionQuantificationView::RA() { TensorQuantify(1); } void QmitkDiffusionQuantificationView::AD() { TensorQuantify(2); } void QmitkDiffusionQuantificationView::RD() { TensorQuantify(3); } void QmitkDiffusionQuantificationView::ClusterAnisotropy() { TensorQuantify(4); } void QmitkDiffusionQuantificationView::MD() { TensorQuantify(5); } void QmitkDiffusionQuantificationView::QBIQuantify(int method) { QBIQuantification(m_QBallImages, method); } void QmitkDiffusionQuantificationView::TensorQuantify(int method) { TensorQuantification(m_TensorImages, method); } void QmitkDiffusionQuantificationView::QBIQuantification( mitk::DataStorage::SetOfObjects::Pointer inImages, int method) { itk::TimeProbe clock; QString status; int nrFiles = inImages->size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { typedef float TOdfPixelType; const int odfsize = QBALL_ODFSIZE; typedef itk::Vector OdfVectorType; typedef itk::Image OdfVectorImgType; mitk::Image* vol = static_cast((*itemiter)->GetData()); OdfVectorImgType::Pointer itkvol = OdfVectorImgType::New(); mitk::CastToItkImage(vol, itkvol); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); float p1 = m_Controls->m_ParamKEdit->text().toFloat(); float p2 = m_Controls->m_ParamPEdit->text().toFloat(); // COMPUTE RA clock.Start(); MBI_INFO << "Computing GFA "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "Computing GFA for %s", nodename.c_str()).toLatin1()); typedef OdfVectorType::ValueType RealValueType; typedef itk::Image< RealValueType, 3 > RAImageType; typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(itkvol); std::string newname; newname.append(nodename); switch(method) { case 0: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); newname.append("GFA"); break; } case 1: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_QUANTILES_HIGH_LOW); newname.append("01"); break; } case 2: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_QUANTILE_HIGH); newname.append("02"); break; } case 3: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_MAX_ODF_VALUE); newname.append("03"); break; } case 4: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_DECONVOLUTION_COEFFS); newname.append("04"); break; } case 5: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_MIN_MAX_NORMALIZED_STANDARD); newname.append("05"); break; } case 6: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_NORMALIZED_ENTROPY); newname.append("06"); break; } case 7: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_NEMATIC_ORDER_PARAMETER); newname.append("07"); break; } case 8: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_QUANTILES_LOW_HIGH); newname.append("08"); break; } case 9: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_QUANTILE_LOW); newname.append("09"); break; } case 10: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_MIN_ODF_VALUE); newname.append("10"); break; } case 11: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_STD_BY_MAX); newname.append("11"); break; } case 12: { p1 = m_Controls->MinAngle->text().toFloat(); p2 = m_Controls->MaxAngle->text().toFloat(); gfaFilter->SetComputationMethod(GfaFilterType::GFA_PRINCIPLE_CURVATURE); QString paramString; paramString = paramString.append("PC%1-%2").arg(p1).arg(p2); newname.append(paramString.toLatin1()); gfaFilter->SetParam1(p1); gfaFilter->SetParam2(p2); break; } case 13: { gfaFilter->SetComputationMethod(GfaFilterType::GFA_GENERALIZED_GFA); QString paramString; paramString = paramString.append("GFAK%1P%2").arg(p1).arg(p2); newname.append(paramString.toLatin1()); gfaFilter->SetParam1(p1); gfaFilter->SetParam2(p2); break; } default: { newname.append("0"); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); } } gfaFilter->Update(); clock.Stop(); typedef itk::Image ImgType; ImgType::Pointer img = ImgType::New(); img->SetSpacing( gfaFilter->GetOutput()->GetSpacing() ); // Set the image spacing img->SetOrigin( gfaFilter->GetOutput()->GetOrigin() ); // Set the image origin img->SetDirection( gfaFilter->GetOutput()->GetDirection() ); // Set the image direction img->SetLargestPossibleRegion( gfaFilter->GetOutput()->GetLargestPossibleRegion()); img->SetBufferedRegion( gfaFilter->GetOutput()->GetLargestPossibleRegion() ); img->Allocate(); itk::ImageRegionIterator ot (img, img->GetLargestPossibleRegion() ); ot.GoToBegin(); itk::ImageRegionConstIterator it (gfaFilter->GetOutput(), gfaFilter->GetOutput()->GetLargestPossibleRegion() ); for (it.GoToBegin(); !it.IsAtEnd(); ++it) { GfaFilterType::OutputImageType::PixelType val = it.Get(); ot.Set(val * m_Controls->m_ScaleImageValuesBox->value()); ++ot; } // GFA TO DATATREE mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk( img.GetPointer() ); image->SetVolume( img->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetProperty( "name", mitk::StringProperty::New(newname) ); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType( mitk::LookupTable::JET_TRANSPARENT ); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable( lut ); node->SetProperty("LookupTable", lut_prop ); GetDataStorage()->Add(node, *itemiter); mitk::StatusBar::GetInstance()->DisplayText("Computation complete."); ++itemiter; } this->GetRenderWindowPart()->RequestUpdate(); } void QmitkDiffusionQuantificationView::TensorQuantification( mitk::DataStorage::SetOfObjects::Pointer inImages, int method) { itk::TimeProbe clock; QString status; int nrFiles = inImages->size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; mitk::Image* vol = static_cast((*itemiter)->GetData()); TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(vol, itkvol); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // COMPUTE FA clock.Start(); MBI_INFO << "Computing FA "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "Computing FA for %s", nodename.c_str()).toLatin1()); typedef itk::Image< TTensorPixelType, 3 > FAImageType; typedef itk::ShiftScaleImageFilter ShiftScaleFilterType; ShiftScaleFilterType::Pointer multi = ShiftScaleFilterType::New(); multi->SetShift(0.0); multi->SetScale(m_Controls->m_ScaleImageValuesBox->value()); typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; if(method == 0) //FA { /* typedef itk::TensorFractionalAnisotropyImageFilter< TensorImageType, FAImageType > FilterType; FilterType::Pointer anisotropyFilter = FilterType::New(); anisotropyFilter->SetInput( itkvol.GetPointer() ); anisotropyFilter->Update(); multi->SetInput(anisotropyFilter->GetOutput()); nodename = QString(nodename.c_str()).append("_FA").toStdString();*/ MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itkvol.GetPointer() ); measurementsCalculator->SetMeasure(MeasurementsType::FA); measurementsCalculator->Update(); multi->SetInput(measurementsCalculator->GetOutput()); nodename = QString(nodename.c_str()).append("_FA").toStdString(); } else if(method == 1) //RA { /*typedef itk::TensorRelativeAnisotropyImageFilter< TensorImageType, FAImageType > FilterType; FilterType::Pointer anisotropyFilter = FilterType::New(); anisotropyFilter->SetInput( itkvol.GetPointer() ); anisotropyFilter->Update(); multi->SetInput(anisotropyFilter->GetOutput()); nodename = QString(nodename.c_str()).append("_RA").toStdString();*/ MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itkvol.GetPointer() ); measurementsCalculator->SetMeasure(MeasurementsType::RA); measurementsCalculator->Update(); multi->SetInput(measurementsCalculator->GetOutput()); nodename = QString(nodename.c_str()).append("_RA").toStdString(); } else if(method == 2) // AD (Axial diffusivity) { MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itkvol.GetPointer() ); measurementsCalculator->SetMeasure(MeasurementsType::AD); measurementsCalculator->Update(); multi->SetInput(measurementsCalculator->GetOutput()); nodename = QString(nodename.c_str()).append("_AD").toStdString(); } else if(method == 3) // RD (Radial diffusivity, (Lambda2+Lambda3)/2 { MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itkvol.GetPointer() ); measurementsCalculator->SetMeasure(MeasurementsType::RD); measurementsCalculator->Update(); multi->SetInput(measurementsCalculator->GetOutput()); nodename = QString(nodename.c_str()).append("_RD").toStdString(); } else if(method == 4) // 1-(Lambda2+Lambda3)/(2*Lambda1) { MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itkvol.GetPointer() ); measurementsCalculator->SetMeasure(MeasurementsType::CA); measurementsCalculator->Update(); multi->SetInput(measurementsCalculator->GetOutput()); nodename = QString(nodename.c_str()).append("_CA").toStdString(); } else if(method == 5) // MD (Mean Diffusivity, (Lambda1+Lambda2+Lambda3)/3 ) { MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); measurementsCalculator->SetInput(itkvol.GetPointer() ); measurementsCalculator->SetMeasure(MeasurementsType::MD); measurementsCalculator->Update(); multi->SetInput(measurementsCalculator->GetOutput()); nodename = QString(nodename.c_str()).append("_MD").toStdString(); } multi->Update(); clock.Stop(); // FA TO DATATREE mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk( multi->GetOutput() ); image->SetVolume( multi->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); node->SetProperty( "name", mitk::StringProperty::New(nodename) ); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); lut->SetType( mitk::LookupTable::JET_TRANSPARENT ); mitk::LookupTableProperty::Pointer lut_prop = mitk::LookupTableProperty::New(); lut_prop->SetLookupTable( lut ); node->SetProperty("LookupTable", lut_prop ); GetDataStorage()->Add(node, *itemiter); ++itemiter; mitk::StatusBar::GetInstance()->DisplayText("Computation complete."); } this->GetRenderWindowPart()->RequestUpdate(); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.h index a705ee0252..96a8432628 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationView.h @@ -1,103 +1,106 @@ /*=================================================================== 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 _QMITKDIFFUSIONQUANTIFICATIONVIEW_H_INCLUDED #define _QMITKDIFFUSIONQUANTIFICATIONVIEW_H_INCLUDED #include #include #include #include "ui_QmitkDiffusionQuantificationViewControls.h" /*! * \ingroup org_mitk_gui_qt_diffusionquantification_internal * * \brief QmitkDiffusionQuantificationView * * Document your class here. */ class QmitkDiffusionQuantificationView : public QmitkAbstractView { friend struct DqSelListener; // this is needed for all Qt objects that should have a MOC object (everything that derives from QObject) Q_OBJECT public: static const std::string VIEW_ID; QmitkDiffusionQuantificationView(); virtual ~QmitkDiffusionQuantificationView(); typedef itk::VectorImage< short, 3 > ItkDwiType; virtual void CreateQtPartControl(QWidget *parent) override; /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; protected slots: void GFACheckboxClicked(); void GFA(); void Curvature(); void FA(); void RA(); void AD(); void RD(); void ClusterAnisotropy(); void MD(); void ADC_DWI(); void MD_DWI(); void QBIQuantify(int method); void QBIQuantification(mitk::DataStorage::SetOfObjects::Pointer inImages, int method) ; void TensorQuantify(int method); void TensorQuantification(mitk::DataStorage::SetOfObjects::Pointer inImages, int method) ; + void DoBallStickCalculation(); + void DoMultiTensorCalculation(); + protected: void DoAdcCalculation(bool fit); /// \brief called by QmitkAbstractView when DataManager's selection has changed virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkDiffusionQuantificationViewControls* m_Controls; mitk::DataStorage::SetOfObjects::Pointer m_QBallImages; mitk::DataStorage::SetOfObjects::Pointer m_TensorImages; mitk::DataStorage::SetOfObjects::Pointer m_DwImages; static const float m_ScaleDAIValues; }; #endif // _QMITKDIFFUSIONQUANTIFICATIONVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationViewControls.ui index 0029b83902..337ab8e6d2 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkDiffusionQuantificationViewControls.ui @@ -1,359 +1,379 @@ QmitkDiffusionQuantificationViewControls 0 0 343 888 0 0 QmitkTemplate Please Select Input Data Image <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true General Parameters Scale Image Values: 9999999.000000000000000 1.000000000000000 Raw diffusion-weighted image false ADC false MD + + + + false + + + Ball-Stick + + + + + + + false + + + Multi-Tensor Fit + + + Q-Ball image QFrame::NoFrame QFrame::Raised 0 0 0 0 Generalized GFA QFrame::NoFrame QFrame::Raised 0 0 0 0 true k true true p true false GFA QFrame::NoFrame QFrame::Raised 0 0 0 0 Min. angle Max. angle false Curvature Tensor image false FA (Fractional Anisotropy) false RA (Relative Anisotropy) false AD (Axial Diffusivity) false RD (Radial Diffusivity) false MD (Mean Diffusivity) false 1-(λ2+λ3)/(2*λ1) Qt::Vertical QSizePolicy::Expanding 20 220 QmitkDataStorageComboBox.h diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp index 20cbffbc41..d8ef803c11 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.cpp @@ -1,1064 +1,1064 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "QmitkTensorReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include #include #include #include #include // itk includes #include "itkTimeProbe.h" //#include "itkTensor.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "mitkTeemDiffusionTensor3DReconstructionImageFilter.h" #include "itkDiffusionTensor3DReconstructionImageFilter.h" #include "itkTensorImageToDiffusionImageFilter.h" #include "itkPointShell.h" #include "itkVector.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkTensorReconstructionWithEigenvalueCorrectionFilter.h" #include "mitkImageCast.h" #include "mitkImageAccessByItk.h" #include #include #include "mitkProperties.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "mitkLookupTableProperty.h" #include "mitkLookupTable.h" #include "mitkImageStatisticsHolder.h" #include #include #include #include #include #include const std::string QmitkTensorReconstructionView::VIEW_ID = "org.mitk.views.tensorreconstruction"; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; using namespace berry; QmitkTensorReconstructionView::QmitkTensorReconstructionView() : QmitkAbstractView(), m_Controls(nullptr) { m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); } QmitkTensorReconstructionView::~QmitkTensorReconstructionView() { } void QmitkTensorReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkTensorReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); Advanced1CheckboxClicked(); } } void QmitkTensorReconstructionView::SetFocus() { m_Controls->m_Advanced1->setFocus(); } void QmitkTensorReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_StartReconstruction), SIGNAL(clicked()), this, SLOT(Reconstruct()) ); connect( (QObject*)(m_Controls->m_Advanced1), SIGNAL(clicked()), this, SLOT(Advanced1CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_TensorsToDWIButton), SIGNAL(clicked()), this, SLOT(TensorsToDWI()) ); connect( (QObject*)(m_Controls->m_TensorsToQbiButton), SIGNAL(clicked()), this, SLOT(TensorsToQbi()) ); connect( (QObject*)(m_Controls->m_ResidualButton), SIGNAL(clicked()), this, SLOT(ResidualCalculation()) ); connect( (QObject*)(m_Controls->m_PerSliceView), SIGNAL(pointSelected(int, int)), this, SLOT(ResidualClicked(int, int)) ); connect( (QObject*)(m_Controls->m_TensorReconstructionThreshold), SIGNAL(valueChanged(int)), this, SLOT(PreviewThreshold(int)) ); } } -void QmitkTensorReconstructionView::ResidualClicked(int slice, itk::SizeValueType volume) +void QmitkTensorReconstructionView::ResidualClicked(int slice, int volume) { // Use image coord to reset crosshair // Find currently selected diffusion image // Update Label // to do: This position should be modified in order to skip B0 volumes that are not taken into account // when calculating residuals // Find the diffusion image mitk::Image* diffImage; mitk::DataNode::Pointer correctNode; mitk::BaseGeometry* geometry; bool isDiffusionImage(false); if (m_DiffusionImage.IsNotNull()) { isDiffusionImage = mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(m_DiffusionImage->GetData())) ; diffImage = static_cast(m_DiffusionImage->GetData()); geometry = m_DiffusionImage->GetData()->GetGeometry(); // Remember the node whose display index must be updated correctNode = mitk::DataNode::New(); correctNode = m_DiffusionImage; } if( isDiffusionImage ) { GradientDirectionContainerType::Pointer dirs = static_cast( diffImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); - for(itk::SizeValueType i=0; iSize() && i<=volume; i++) + for(itk::SizeValueType i=0; iSize() && i<=(itk::SizeValueType)volume; i++) { GradientDirectionType grad = dirs->ElementAt(i); // check if image is b0 weighted if(fabs(grad[0]) < 0.001 && fabs(grad[1]) < 0.001 && fabs(grad[2]) < 0.001) { volume++; } } QString pos = "Volume: "; pos.append(QString::number(volume)); pos.append(", Slice: "); pos.append(QString::number(slice)); m_Controls->m_PositionLabel->setText(pos); if(correctNode) { int oldDisplayVal; correctNode->GetIntProperty("DisplayChannel", oldDisplayVal); QString oldVal = QString::number(oldDisplayVal); QString newVal = QString::number(volume); correctNode->SetIntProperty("DisplayChannel",volume); correctNode->SetSelected(true); this->FirePropertyChanged("DisplayChannel", oldVal, newVal); correctNode->UpdateOutputInformation(); mitk::Point3D p3 = this->GetRenderWindowPart()->GetSelectedPosition(); itk::Index<3> ix; geometry->WorldToIndex(p3, ix); // ix[2] = slice; mitk::Vector3D vec; vec[0] = ix[0]; vec[1] = ix[1]; vec[2] = slice; mitk::Vector3D v3New; geometry->IndexToWorld(vec, v3New); mitk::Point3D origin = geometry->GetOrigin(); mitk::Point3D p3New; p3New[0] = v3New[0] + origin[0]; p3New[1] = v3New[1] + origin[1]; p3New[2] = v3New[2] + origin[2]; this->GetRenderWindowPart()->SetSelectedPosition(p3New); this->GetRenderWindowPart()->RequestUpdate(); } } } void QmitkTensorReconstructionView::Advanced1CheckboxClicked() { bool check = m_Controls-> m_Advanced1->isChecked(); m_Controls->frame->setVisible(check); } void QmitkTensorReconstructionView::Activated() { } void QmitkTensorReconstructionView::Deactivated() { // Get all current nodes mitk::DataStorage::SetOfObjects::ConstPointer objects = this->GetDataStorage()->GetAll(); mitk::DataStorage::SetOfObjects::const_iterator itemiter( objects->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( objects->end() ); while ( itemiter != itemiterend ) // for all items { mitk::DataNode::Pointer node = *itemiter; if (node.IsNull()) continue; // only look at interesting types if( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage(dynamic_cast(node->GetData()))) { if (this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", *itemiter)) { node = this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", *itemiter); this->GetDataStorage()->Remove(node); } } itemiter++; } mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkTensorReconstructionView::Visible() { } void QmitkTensorReconstructionView::Hidden() { } void QmitkTensorReconstructionView::ResidualCalculation() { // Extract dwi and dti from current selection // In case of multiple selections, take the first one, since taking all combinations is not meaningful mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); mitk::Image::Pointer diffImage = mitk::Image::New(); TensorImageType::Pointer tensorImage; std::string nodename; if(m_DiffusionImage.IsNotNull()) { diffImage = static_cast(m_DiffusionImage->GetData()); } else return; if(m_TensorImage.IsNotNull()) { mitk::TensorImage* mitkVol; mitkVol = static_cast(m_TensorImage->GetData()); mitk::CastToItkImage(mitkVol, tensorImage); m_TensorImage->GetStringProperty("name", nodename); } else return; typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; GradientDirectionContainerType* gradients = static_cast( diffImage->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer(); // Find the min and the max values from a baseline image mitk::ImageStatisticsHolder *stats = diffImage->GetStatistics(); //Initialize filter that calculates the modeled diffusion weighted signals FilterType::Pointer filter = FilterType::New(); filter->SetInput( tensorImage ); filter->SetBValue( static_cast(diffImage->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue()); filter->SetGradientList(gradients); filter->SetMin(stats->GetScalarValueMin()); filter->SetMax(stats->GetScalarValueMax()); filter->Update(); // TENSORS TO DATATREE mitk::Image::Pointer image = mitk::GrabItkImageMemory( filter->GetOutput() ); image->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( gradients ) ); image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( static_cast(diffImage->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ) ); image->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New( static_cast(diffImage->GetProperty(mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str()).GetPointer() )->GetMeasurementFrame() ) ); mitk::DiffusionPropertyHelper propertyHelper( image ); propertyHelper.InitializeImage(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); mitk::ImageVtkMapper2D::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_Estimated DWI"); node->SetName(newname.toLatin1()); GetDataStorage()->Add(node, m_TensorImage); BValueMapType map = static_cast(image->GetProperty(mitk::DiffusionPropertyHelper::BVALUEMAPPROPERTYNAME.c_str()).GetPointer() )->GetBValueMap();; std::vector< unsigned int > b0Indices = map[0]; typedef itk::ResidualImageFilter ResidualImageFilterType; ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(diffImage, itkVectorImagePointer); ITKDiffusionImageType::Pointer itkSecondVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(image, itkSecondVectorImagePointer); ResidualImageFilterType::Pointer residualFilter = ResidualImageFilterType::New(); residualFilter->SetInput( itkVectorImagePointer ); residualFilter->SetSecondDiffusionImage( itkSecondVectorImagePointer ); residualFilter->SetGradients(gradients); residualFilter->SetB0Index(b0Indices[0]); residualFilter->SetB0Threshold(30); residualFilter->Update(); itk::Image::Pointer residualImage = itk::Image::New(); residualImage = residualFilter->GetOutput(); mitk::Image::Pointer mitkResImg = mitk::Image::New(); mitk::CastToMitkImage(residualImage, mitkResImg); stats = mitkResImg->GetStatistics(); float min = stats->GetScalarValueMin(); float max = stats->GetScalarValueMax(); mitk::LookupTableProperty::Pointer lutProp = mitk::LookupTableProperty::New(); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(min, max); // If you don't want to use the whole color range, you can use // SetValueRange, SetHueRange, and SetSaturationRange lookupTable->Build(); vtkSmartPointer reversedlookupTable = vtkSmartPointer::New(); reversedlookupTable->SetTableRange(min+1, max); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookupTable->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } lut->SetVtkLookupTable(lookupTable); lutProp->SetLookupTable(lut); // Create lookuptable mitk::DataNode::Pointer resNode=mitk::DataNode::New(); resNode->SetData( mitkResImg ); resNode->SetName("Residuals"); resNode->SetProperty("LookupTable", lutProp); bool b; resNode->GetBoolProperty("use color", b); resNode->SetBoolProperty("use color", false); GetDataStorage()->Add(resNode, m_TensorImage); this->GetRenderWindowPart()->RequestUpdate(); // Draw Graph std::vector means = residualFilter->GetMeans(); std::vector q1s = residualFilter->GetQ1(); std::vector q3s = residualFilter->GetQ3(); std::vector percentagesOfOUtliers = residualFilter->GetPercentagesOfOutliers(); m_Controls->m_ResidualAnalysis->SetMeans(means); m_Controls->m_ResidualAnalysis->SetQ1(q1s); m_Controls->m_ResidualAnalysis->SetQ3(q3s); m_Controls->m_ResidualAnalysis->SetPercentagesOfOutliers(percentagesOfOUtliers); if(m_Controls->m_PercentagesOfOutliers->isChecked()) { m_Controls->m_ResidualAnalysis->DrawPercentagesOfOutliers(); } else { m_Controls->m_ResidualAnalysis->DrawMeans(); } // Draw Graph for volumes per slice in the QGraphicsView std::vector< std::vector > outliersPerSlice = residualFilter->GetOutliersPerSlice(); int xSize = outliersPerSlice.size(); if(xSize == 0) { return; } int ySize = outliersPerSlice[0].size(); // Find maximum in outliersPerSlice double maxOutlier= 0.0; for(int i=0; imaxOutlier) { maxOutlier = outliersPerSlice[i][j]; } } } // Create some QImage QImage qImage(xSize, ySize, QImage::Format_RGB32); QImage legend(1, 256, QImage::Format_RGB32); QRgb value; vtkSmartPointer lookup = vtkSmartPointer::New(); lookup->SetTableRange(0.0, maxOutlier); lookup->Build(); reversedlookupTable->SetTableRange(0, maxOutlier); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookup->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } // Fill qImage for(int i=0; iMapValue(out); int r, g, b; r = _rgba[0]; g = _rgba[1]; b = _rgba[2]; value = qRgb(r, g, b); qImage.setPixel(i,j,value); } } for(int i=0; i<256; i++) { double* rgba = lookup->GetTableValue(i); int r, g, b; r = rgba[0]*255; g = rgba[1]*255; b = rgba[2]*255; value = qRgb(r, g, b); legend.setPixel(0,255-i,value); } QString upper = QString::number(maxOutlier, 'g', 3); upper.append(" %"); QString lower = QString::number(0.0); lower.append(" %"); m_Controls->m_UpperLabel->setText(upper); m_Controls->m_LowerLabel->setText(lower); QPixmap pixmap(QPixmap::fromImage(qImage)); QGraphicsPixmapItem *item = new QGraphicsPixmapItem(pixmap); item->setTransform(QTransform::fromScale(10.0, 3.0), true); QPixmap pixmap2(QPixmap::fromImage(legend)); QGraphicsPixmapItem *item2 = new QGraphicsPixmapItem(pixmap2); item2->setTransform(QTransform::fromScale(20.0, 1.0), true); m_Controls->m_PerSliceView->SetResidualPixmapItem(item); QGraphicsScene* scene = new QGraphicsScene; QGraphicsScene* scene2 = new QGraphicsScene; scene->addItem(item); scene2->addItem(item2); m_Controls->m_PerSliceView->setScene(scene); m_Controls->m_LegendView->setScene(scene2); m_Controls->m_PerSliceView->show(); m_Controls->m_PerSliceView->repaint(); m_Controls->m_LegendView->setHorizontalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->setVerticalScrollBarPolicy ( Qt::ScrollBarAlwaysOff ); m_Controls->m_LegendView->show(); m_Controls->m_LegendView->repaint(); } void QmitkTensorReconstructionView::Reconstruct() { int method = m_Controls->m_ReconctructionMethodBox->currentIndex(); switch (method) { case 0: ItkTensorReconstruction(m_DiffusionImages); break; case 1: TensorReconstructionWithCorr(m_DiffusionImages); break; default: ItkTensorReconstruction(m_DiffusionImages); } } void QmitkTensorReconstructionView::TensorReconstructionWithCorr (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); while ( itemiter != itemiterend ) // for all items { mitk::Image* vols = static_cast((*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // TENSOR RECONSTRUCTION MITK_INFO << "Tensor reconstruction with correction for negative eigenvalues"; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Tensor reconstruction for %s", nodename.c_str()).toLatin1()); typedef itk::TensorReconstructionWithEigenvalueCorrectionFilter< DiffusionPixelType, TTensorPixelType > ReconstructionFilter; float b0Threshold = m_Controls->m_TensorReconstructionThreshold->value(); GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = static_cast( vols->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->Begin(); it != static_cast( vols->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->End(); it++) { gradientContainerCopy->push_back(it.Value()); } ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); ReconstructionFilter::Pointer reconFilter = ReconstructionFilter::New(); reconFilter->SetBValue( static_cast(vols->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); reconFilter->SetGradientImage( gradientContainerCopy, itkVectorImagePointer); reconFilter->SetB0Threshold(b0Threshold); reconFilter->Update(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer outputTensorImg = reconFilter->GetOutput(); typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(outputTensorImg, outputTensorImg->GetRequestedRegion()); tensorIt.GoToBegin(); int negatives = 0; while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iInitializeByItk( outputTensorImg.GetPointer() ); image->SetVolume( outputTensorImg->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); SetDefaultNodeProperties(node, nodename+"_EigenvalueCorrected_DT"); GetDataStorage()->Add(node, *itemiter); mitk::ProgressBar::GetInstance()->Progress(); ++itemiter; } mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toLatin1()); this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); } } void QmitkTensorReconstructionView::ItkTensorReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); while ( itemiter != itemiterend ) // for all items { mitk::Image* vols = static_cast( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // TENSOR RECONSTRUCTION clock.Start(); MITK_DEBUG << "Tensor reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Tensor reconstruction for %s", nodename.c_str()).toLatin1()); typedef itk::DiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); GradientDirectionContainerType::Pointer gradientContainerCopy = GradientDirectionContainerType::New(); for(GradientDirectionContainerType::ConstIterator it = static_cast( vols->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->Begin(); it != static_cast( vols->GetProperty(mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str()).GetPointer() )->GetGradientDirectionsContainer()->End(); it++) { gradientContainerCopy->push_back(it.Value()); } ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); tensorReconstructionFilter->SetBValue( static_cast(vols->GetProperty(mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str()).GetPointer() )->GetValue() ); tensorReconstructionFilter->SetGradientImage( gradientContainerCopy, itkVectorImagePointer ); tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreshold->value() ); tensorReconstructionFilter->Update(); clock.Stop(); MITK_DEBUG << "took " << clock.GetMean() << "s."; // TENSORS TO DATATREE mitk::TensorImage::Pointer image = mitk::TensorImage::New(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer tensorImage; tensorImage = tensorReconstructionFilter->GetOutput(); // Check the tensor for negative eigenvalues if(m_Controls->m_CheckNegativeEigenvalues->isChecked()) { typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetRequestedRegion()); tensorIt.GoToBegin(); while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; //typedef itk::Tensor TensorType2; TensorType tensor = tensorIt.Get(); TensorType::EigenValuesArrayType ev; tensor.ComputeEigenValues(ev); for(unsigned int i=0; iSetDirection( itkVectorImagePointer->GetDirection() ); image->InitializeByItk( tensorImage.GetPointer() ); image->SetVolume( tensorReconstructionFilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); SetDefaultNodeProperties(node, nodename+"_LinearLeastSquares_DT"); GetDataStorage()->Add(node, *itemiter); mitk::ProgressBar::GetInstance()->Progress(); ++itemiter; } mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toLatin1()); this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return; } } void QmitkTensorReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name) { node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) ); node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) ); node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New()); node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New()); node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2)); node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1)); node->SetProperty( "visible", mitk::BoolProperty::New( true ) ); node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) ); node->SetProperty ("layer", mitk::IntProperty::New(100)); node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) ); node->SetProperty( "name", mitk::StringProperty::New(name) ); } void QmitkTensorReconstructionView::TensorsToDWI() { DoTensorsToDWI(m_TensorImages); } void QmitkTensorReconstructionView::TensorsToQbi() { for (unsigned int i=0; isize(); i++) { mitk::DataNode::Pointer tensorImageNode = m_TensorImages->at(i); MITK_INFO << "starting Q-Ball estimation"; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(dynamic_cast(tensorImageNode->GetData()), itkvol); typedef itk::TensorImageToQBallImageFilter< TTensorPixelType, TTensorPixelType > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->Update(); typedef itk::Vector OutputPixelType; typedef itk::Image OutputImageType; mitk::QBallImage::Pointer image = mitk::QBallImage::New(); OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); node->SetName(tensorImageNode->GetName()+"_Qball"); GetDataStorage()->Add(node, tensorImageNode); } } void QmitkTensorReconstructionView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { m_DiffusionImages = mitk::DataStorage::SetOfObjects::New(); m_TensorImages = mitk::DataStorage::SetOfObjects::New(); bool foundDwiVolume = false; bool foundTensorVolume = false; m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_DiffusionImage = nullptr; m_TensorImage = nullptr; m_Controls->m_InputData->setTitle("Please Select Input Data"); // iterate selection for (mitk::DataNode::Pointer node: nodes) { if (node.IsNull()) continue; // only look at interesting types bool isDiffusionImage( mitk::DiffusionPropertyHelper::IsDiffusionWeightedImage( dynamic_cast(node->GetData())) ); if ( isDiffusionImage ) { foundDwiVolume = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_DiffusionImages->push_back(node); m_DiffusionImage = node; } else if(dynamic_cast(node->GetData())) { foundTensorVolume = true; m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_TensorImages->push_back(node); m_TensorImage = node; } } m_Controls->m_StartReconstruction->setEnabled(foundDwiVolume); m_Controls->m_TensorsToDWIButton->setEnabled(foundTensorVolume); m_Controls->m_TensorsToQbiButton->setEnabled(foundTensorVolume); if (foundDwiVolume || foundTensorVolume) m_Controls->m_InputData->setTitle("Input Data"); m_Controls->m_ResidualButton->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PercentagesOfOutliers->setEnabled(foundDwiVolume && foundTensorVolume); m_Controls->m_PerSliceView->setEnabled(foundDwiVolume && foundTensorVolume); } template itk::VectorContainer >::Pointer QmitkTensorReconstructionView::MakeGradientList() { itk::VectorContainer >::Pointer retval = itk::VectorContainer >::New(); vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); for(int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval->push_back(v); } // Add 0 vector for B0 vnl_vector_fixed v; v.fill(0.0); retval->push_back(v); return retval; } void QmitkTensorReconstructionView::DoTensorsToDWI(mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); while ( itemiter != itemiterend ) // for all items { std::string nodename; (*itemiter)->GetStringProperty("name", nodename); mitk::TensorImage* vol = static_cast((*itemiter)->GetData()); typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(vol, itkvol); typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; FilterType::GradientListPointerType gradientList = FilterType::GradientListType::New(); switch(m_Controls->m_TensorsToDWINumDirsSelect->currentIndex()) { case 0: gradientList = MakeGradientList<12>(); break; case 1: gradientList = MakeGradientList<42>(); break; case 2: gradientList = MakeGradientList<92>(); break; case 3: gradientList = MakeGradientList<162>(); break; case 4: gradientList = MakeGradientList<252>(); break; case 5: gradientList = MakeGradientList<362>(); break; case 6: gradientList = MakeGradientList<492>(); break; case 7: gradientList = MakeGradientList<642>(); break; case 8: gradientList = MakeGradientList<812>(); break; case 9: gradientList = MakeGradientList<1002>(); break; default: gradientList = MakeGradientList<92>(); } double bVal = m_Controls->m_TensorsToDWIBValueEdit->text().toDouble(); // DWI ESTIMATION clock.Start(); MBI_INFO << "DWI Estimation "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "DWI Estimation for %s", nodename.c_str()).toLatin1()); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->SetBValue(bVal); filter->SetGradientList(gradientList); //filter->SetNumberOfThreads(1); filter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMean() << "s."; // TENSORS TO DATATREE mitk::Image::Pointer image = mitk::GrabItkImageMemory( filter->GetOutput() ); image->SetProperty( mitk::DiffusionPropertyHelper::GRADIENTCONTAINERPROPERTYNAME.c_str(), mitk::GradientDirectionsProperty::New( gradientList ) ); image->SetProperty( mitk::DiffusionPropertyHelper::REFERENCEBVALUEPROPERTYNAME.c_str(), mitk::FloatProperty::New( bVal ) ); image->SetProperty( mitk::DiffusionPropertyHelper::MEASUREMENTFRAMEPROPERTYNAME.c_str(), mitk::MeasurementFrameProperty::New() ); mitk::DiffusionPropertyHelper propertyHelper( image ); propertyHelper.InitializeImage(); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); mitk::ImageVtkMapper2D::SetDefaultProperties(node); node->SetName(nodename+"_DWI"); GetDataStorage()->Add(node, *itemiter); mitk::ProgressBar::GetInstance()->Progress(); ++itemiter; } mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toLatin1()); this->GetRenderWindowPart()->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "DWI estimation failed:", ex.GetDescription()); return ; } } void QmitkTensorReconstructionView::PreviewThreshold(int threshold) { mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_DiffusionImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_DiffusionImages->end() ); while ( itemiter != itemiterend ) // for all items { mitk::Image* vols = static_cast( (*itemiter)->GetData()); ITKDiffusionImageType::Pointer itkVectorImagePointer = ITKDiffusionImageType::New(); mitk::CastToItkImage(vols, itkVectorImagePointer); // Extract b0 image typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filterB0 = FilterType::New(); filterB0->SetInput(itkVectorImagePointer); filterB0->SetDirections(mitk::DiffusionPropertyHelper::GetGradientContainer(vols)); filterB0->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); typedef itk::Image ImageType; typedef itk::Image SegmentationType; typedef itk::BinaryThresholdImageFilter ThresholdFilterType; // apply threshold ThresholdFilterType::Pointer filterThreshold = ThresholdFilterType::New(); filterThreshold->SetInput(filterB0->GetOutput()); filterThreshold->SetLowerThreshold(threshold); filterThreshold->SetInsideValue(0); filterThreshold->SetOutsideValue(1); // mark cut off values red filterThreshold->Update(); mitkImage->InitializeByItk( filterThreshold->GetOutput() ); mitkImage->SetVolume( filterThreshold->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node; if (this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", *itemiter)) { node = this->GetDataStorage()->GetNamedDerivedNode("ThresholdOverlay", *itemiter); } else { // create a new node, to show thresholded values node = mitk::DataNode::New(); GetDataStorage()->Add( node, *itemiter ); node->SetProperty( "name", mitk::StringProperty::New("ThresholdOverlay")); node->SetBoolProperty("helper object", true); } node->SetData( mitkImage ); itemiter++; mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.h b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.h index 626349b333..eb83870ab5 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.reconstruction/src/internal/QmitkTensorReconstructionView.h @@ -1,130 +1,130 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _QMITKTENSORRECONSTRUCTIONVIEW_H_INCLUDED #define _QMITKTENSORRECONSTRUCTIONVIEW_H_INCLUDED #include #include #include #include "ui_QmitkTensorReconstructionViewControls.h" #include #include #include #include typedef short DiffusionPixelType; struct TrSelListener; /*! * \ingroup org_mitk_gui_qt_tensorreconstruction_internal * * \brief QmitkTensorReconstructionView * * Document your class here. */ class QmitkTensorReconstructionView : public QmitkAbstractView, public mitk::ILifecycleAwarePart { friend struct TrSelListener; // this is needed for all Qt objects that should have a MOC object (everything that derives from QObject) Q_OBJECT public: typedef mitk::DiffusionPropertyHelper::GradientDirectionType GradientDirectionType; typedef mitk::DiffusionPropertyHelper::GradientDirectionsContainerType GradientDirectionContainerType; typedef mitk::DiffusionPropertyHelper::BValueMapType BValueMapType; typedef itk::VectorImage< DiffusionPixelType, 3 > ITKDiffusionImageType; static const std::string VIEW_ID; QmitkTensorReconstructionView(); virtual ~QmitkTensorReconstructionView(); virtual void CreateQtPartControl(QWidget *parent) override; /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); /// /// Sets the focus to an internal widget. /// virtual void SetFocus() override; /// \brief Called when the view gets activated virtual void Activated() override; /// \brief Called when the view gets deactivated virtual void Deactivated() override; /// \brief Called when the view becomes visible virtual void Visible() override; /// \brief Called when the view becomes hidden virtual void Hidden() override; static const int nrconvkernels; protected slots: void TensorsToQbi(); void TensorsToDWI(); void DoTensorsToDWI(mitk::DataStorage::SetOfObjects::Pointer inImages); void Advanced1CheckboxClicked(); void Reconstruct(); void ResidualCalculation(); - void ResidualClicked(int slice, itk::SizeValueType volume); + void ResidualClicked(int slice, int volume); /** * @brief PreviewThreshold Generates a preview of the values that are cut off by the thresholds * @param threshold */ void PreviewThreshold(int threshold); protected: void ItkTensorReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages); void TeemTensorReconstruction(mitk::DataStorage::SetOfObjects::Pointer inImages); void TensorReconstructionWithCorr(mitk::DataStorage::SetOfObjects::Pointer inImages); virtual void OnSelectionChanged(berry::IWorkbenchPart::Pointer part, const QList& nodes) override; Ui::QmitkTensorReconstructionViewControls* m_Controls; template itk::VectorContainer >::Pointer MakeGradientList(); template void TemplatedAnalyticalTensorReconstruction(mitk::Image* vols, float lambda, std::string nodename, std::vector* nodes, int normalization); void SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name); mitk::DataNode::Pointer m_DiffusionImage; mitk::DataNode::Pointer m_TensorImage; mitk::DataStorage::SetOfObjects::Pointer m_DiffusionImages; mitk::DataStorage::SetOfObjects::Pointer m_TensorImages; }; #endif // _QMITKTENSORRECONSTRUCTIONVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/QmitkDiffusionImagingAppWorkbenchAdvisor.cpp b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/QmitkDiffusionImagingAppWorkbenchAdvisor.cpp index e72a83ad70..26dd08c6ca 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/QmitkDiffusionImagingAppWorkbenchAdvisor.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimagingapp/src/QmitkDiffusionImagingAppWorkbenchAdvisor.cpp @@ -1,99 +1,99 @@ /*=================================================================== 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 "QmitkDiffusionImagingAppWorkbenchAdvisor.h" #include "internal/QmitkDiffusionApplicationPlugin.h" #include #include #include #include #include #include #include #include const QString QmitkDiffusionImagingAppWorkbenchAdvisor::WELCOME_PERSPECTIVE_ID = "org.mitk.diffusionimagingapp.perspectives.welcome"; class QmitkDiffusionImagingAppWorkbenchWindowAdvisor : public QmitkExtWorkbenchWindowAdvisor { public: QmitkDiffusionImagingAppWorkbenchWindowAdvisor(berry::WorkbenchAdvisor* wbAdvisor, berry::IWorkbenchWindowConfigurer::Pointer configurer) : QmitkExtWorkbenchWindowAdvisor(wbAdvisor, configurer) { } void PostWindowOpen() { QmitkExtWorkbenchWindowAdvisor::PostWindowOpen(); berry::IWorkbenchWindowConfigurer::Pointer configurer = GetWindowConfigurer(); configurer->GetWindow()->GetWorkbench()->GetIntroManager()->ShowIntro(configurer->GetWindow(), false); } }; void QmitkDiffusionImagingAppWorkbenchAdvisor::Initialize(berry::IWorkbenchConfigurer::Pointer configurer) { berry::QtWorkbenchAdvisor::Initialize(configurer); configurer->SetSaveAndRestore(true); } berry::WorkbenchWindowAdvisor* QmitkDiffusionImagingAppWorkbenchAdvisor::CreateWorkbenchWindowAdvisor( berry::IWorkbenchWindowConfigurer::Pointer configurer) { QList perspExcludeList; perspExcludeList.push_back( "org.blueberry.uitest.util.EmptyPerspective" ); perspExcludeList.push_back( "org.blueberry.uitest.util.EmptyPerspective2" ); perspExcludeList.push_back("org.blueberry.perspectives.help"); QList viewExcludeList; viewExcludeList.push_back( "org.mitk.views.controlvisualizationpropertiesview" ); viewExcludeList.push_back( "org.mitk.views.modules" ); - QRect rec = QApplication::desktop()->screenGeometry(); - configurer->SetInitialSize(QPoint(rec.width(),rec.height())); + //QRect rec = QApplication::desktop()->screenGeometry(); + //configurer->SetInitialSize(QPoint(rec.width(),rec.height())); QmitkDiffusionImagingAppWorkbenchWindowAdvisor* advisor = new QmitkDiffusionImagingAppWorkbenchWindowAdvisor(this, configurer); advisor->ShowViewMenuItem(true); advisor->ShowNewWindowMenuItem(true); advisor->ShowClosePerspectiveMenuItem(true); advisor->SetPerspectiveExcludeList(perspExcludeList); advisor->SetViewExcludeList(viewExcludeList); advisor->ShowViewToolbar(false); advisor->ShowPerspectiveToolbar(true); advisor->ShowVersionInfo(false); advisor->ShowMitkVersionInfo(false); advisor->ShowMemoryIndicator(false); advisor->SetProductName("MITK Diffusion"); advisor->SetWindowIcon(":/org.mitk.gui.qt.diffusionimagingapp/app-icon.png"); return advisor; } QString QmitkDiffusionImagingAppWorkbenchAdvisor::GetInitialWindowPerspectiveId() { return WELCOME_PERSPECTIVE_ID; }