diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h index c991531c29..d191ab7945 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h @@ -1,219 +1,219 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkDiffusionMultiShellQballReconstructionImageFilter_h_ #define __itkDiffusionMultiShellQballReconstructionImageFilter_h_ #include namespace itk{ /** \class DiffusionMultiShellQballReconstructionImageFilter I. Aganj, C. Lenglet, G. Sapiro, E. Yacoub, K. Ugurbil, and N. Harel, “Reconstruction of the orientation distribution function in single and multiple shell q-ball imaging within constant solid angle,” Magnetic Resonance in Medicine, vol. 64, no. 2, pp. 554–566, 2010. */ template< class TReferenceImagePixelType, class TGradientImagePixelType, class TOdfPixelType, int NOrderL, int NrOdfDirections> class DiffusionMultiShellQballReconstructionImageFilter : public ImageToImageFilter< Image< TReferenceImagePixelType, 3 >, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > { public: typedef DiffusionMultiShellQballReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< TReferenceImagePixelType, 3>, Image< Vector< TOdfPixelType, NrOdfDirections >, 3 > > Superclass; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef TReferenceImagePixelType ReferencePixelType; /** GradientImageType * (e.g. type short)*/ typedef TGradientImagePixelType GradientPixelType; /** GradientImageType * 3D VectorImage containing GradientPixelTypes */ typedef VectorImage< GradientPixelType, 3 > GradientImagesType; /** ODF PixelType */ typedef Vector< TOdfPixelType, NrOdfDirections > OdfPixelType; /** ODF ImageType */ typedef Image< OdfPixelType, 3 > OdfImageType; /** BzeroImageType */ typedef Image< TOdfPixelType, 3 > BZeroImageType; /** Container to hold gradient directions of the 'n' DW measurements */ typedef VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType; typedef Image< Vector< TOdfPixelType, (NOrderL*NOrderL + NOrderL + 2)/2 + NOrderL >, 3 > CoefficientImageType; - typedef std::map > BValueMap; - typedef std::map >::iterator BValueMapIteraotr; + typedef std::map > BValueMap; + typedef std::map >::iterator BValueMapIteraotr; typedef std::vector IndiciesVector; // --------------------------------------------------------------------------------------------// /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(DiffusionMultiShellQballReconstructionImageFilter, ImageToImageFilter) /** Get reference image */ virtual typename Superclass::InputImageType * GetInputImage() { return ( static_cast< typename Superclass::InputImageType *>(this->ProcessObject::GetInput(0)) ); } /** Replaces the Input method. * Var vols = mitk::DiffusionImage * ----------------------------------------------------- * GradientDirectionContainerType-Input gradientDirectionContainer (e.g. vols->GetDirections) * GradientImagesType-Input gradientImage (e.g. vols->GetVectorImage) * float-Input bvalue (e.g. vols->GetB_Value) */ void SetGradientImage( const GradientDirectionContainerType * gradientDirectionContainer, const GradientImagesType *gradientImage , float bvalue);//, std::vector listOfUserSelctedBValues ); /** Set a BValue Map (key = bvalue, value = indicies splittet for each shell) * If the input image containes more than three q-shells * (e.g. b-Values of 0, 1000, 2000, 3000, 4000, ...). * For the Analytical-Reconstruction it is needed to set a * BValue Map containing three shells in an arithmetic series * (e.g. 0, 1000, 2000, 3000). */ inline void SetBValueMap(BValueMap map){this->m_BValueMap = map;} /** Threshold on the reference image data. The output ODF will be a null * pdf for pixels in the reference image that have a value less than this * threshold. */ itkSetMacro( Threshold, ReferencePixelType ) itkGetMacro( Threshold, ReferencePixelType ) itkGetMacro( CoefficientImage, typename CoefficientImageType::Pointer ) /** Return non-diffusion weighted images */ itkGetMacro( BZeroImage, typename BZeroImageType::Pointer) /** Factor for Laplacian-Baltrami smoothing of the SH-coefficients*/ itkSetMacro( Lambda, double ) itkGetMacro( Lambda, double ) protected: DiffusionMultiShellQballReconstructionImageFilter(); ~DiffusionMultiShellQballReconstructionImageFilter() { } void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads ); private: enum ReconstructionType { Mode_Analytical3Shells, Mode_NumericalNShells, Mode_Standard1Shell }; ReconstructionType m_ReconstructionType; // Interpolation bool m_Interpolation_Flag; vnl_matrix< double > * m_Interpolation_SHT1_inv; vnl_matrix< double > * m_Interpolation_SHT2_inv; vnl_matrix< double > * m_Interpolation_SHT3_inv; vnl_matrix< double > * m_TARGET_SH_shell1; vnl_matrix< double > * m_TARGET_SH_shell2; vnl_matrix< double > * m_TARGET_SH_shell3; unsigned int m_MaxDirections; vnl_matrix< double > * m_CoeffReconstructionMatrix; vnl_matrix< double > * m_ODFSphericalHarmonicBasisMatrix; /** container to hold gradient directions */ GradientDirectionContainerType::Pointer m_GradientDirectionContainer; /** Number of gradient measurements */ unsigned int m_NumberOfGradientDirections; /** Number of baseline images */ unsigned int m_NumberOfBaselineImages; /** Threshold on the reference image data */ ReferencePixelType m_Threshold; typename BZeroImageType::Pointer m_BZeroImage; typename CoefficientImageType::Pointer m_CoefficientImage; float m_BValue; BValueMap m_BValueMap; double m_Lambda; bool m_IsHemisphericalArrangementOfGradientDirections; bool m_IsArithmeticProgession; void ComputeReconstructionMatrix(IndiciesVector const & refVector); void ComputeODFSHBasis(); bool CheckDuplicateDiffusionGradients(); bool CheckForDifferingShellDirections(); IndiciesVector GetAllDirections(); void ComputeSphericalHarmonicsBasis(vnl_matrix* QBallReference, vnl_matrix* SHBasisOutput, int Lorder , vnl_matrix* LaplaciaBaltramiOutput =0 , vnl_vector* SHOrderAssociation =0 , vnl_matrix * SHEigenvalues =0); void Normalize(OdfPixelType & odf ); void S_S0Normalization( vnl_vector & vec, double b0 = 0 ); void DoubleLogarithm(vnl_vector & vec); double CalculateThreashold(const double value, const double delta); void Projection1(vnl_vector & vec, double delta = 0.01); void Projection2( vnl_vector & E1, vnl_vector & E2, vnl_vector & E3, double delta = 0.01); void Projection3( vnl_vector & A, vnl_vector & alpha, vnl_vector & beta, double delta = 0.01); void StandardOneShellReconstruction(const OutputImageRegionType& outputRegionForThread); void AnalyticalThreeShellReconstruction(const OutputImageRegionType& outputRegionForThread); void NumericalNShellReconstruction(const OutputImageRegionType& outputRegionForThread); void GenerateAveragedBZeroImage(const OutputImageRegionType& outputRegionForThread); void ComputeSphericalFromCartesian(vnl_matrix * Q, const IndiciesVector & refShell); //------------------------- VNL-function ------------------------------------ template vnl_vector< WntValue> element_cast (vnl_vector< CurrentValue> const& v1) { vnl_vector result(v1.size()); for(unsigned int i = 0 ; i < v1.size(); i++) result[i] = static_cast< WntValue>(v1[i]); return result; } template double dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ) { double result = (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (v1.two_norm() * v2.two_norm()); return result ; } }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDiffusionMultiShellQballReconstructionImageFilter.cpp" #endif #endif //__itkDiffusionMultiShellQballReconstructionImageFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h index fc641277ef..5af09ee771 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h @@ -1,119 +1,119 @@ /*=================================================================== 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 _itk_MultiShellAdcAverageReconstructionImageFilter_h_ #define _itk_MultiShellAdcAverageReconstructionImageFilter_h_ #include #include #include namespace itk { /** * \brief Select subset of the input vectors equally distributed over the sphere using an iterative electrostatic repulsion strategy. */ template class MultiShellAdcAverageReconstructionImageFilter : public ImageToImageFilter, itk::VectorImage > { public: typedef MultiShellAdcAverageReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< itk::VectorImage, itk::VectorImage > Superclass; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(MultiShellAdcAverageReconstructionImageFilter, ImageToImageFilter) typedef TInputScalarType InputScalarType; typedef itk::VectorImage InputImageType; typedef typename InputImageType::PixelType InputPixelType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef std::vector IndicesVector; - typedef std::map BValueMap; + typedef std::map BValueMap; GradientDirectionContainerType::Pointer GetOriginalGradientDirections(){return m_OriginalGradientDirections;} void SetOriginalGradientDirections(GradientDirectionContainerType::Pointer ptr){m_OriginalGradientDirections = ptr;} GradientDirectionContainerType::Pointer GetTargetGradientDirections(){return m_TargetGradientDirections;} double GetTargetB_Value(){return m_TargetB_Value;} double GetB_Value(){return m_B_Value;} void SetB_Value(double val){m_B_Value = val;} void SetOriginalBValueMap(BValueMap inp){m_B_ValueMap = inp;} protected: MultiShellAdcAverageReconstructionImageFilter(); ~MultiShellAdcAverageReconstructionImageFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads ); void S_S0Normalization( vnl_vector & vec, const double & S0 ); void calculateAdcFromSignal( vnl_vector & vec, const double & bValue); void calculateSignalFromAdc( vnl_vector & vec, const double & bValue, const double & referenceSignal); GradientDirectionContainerType::Pointer m_TargetGradientDirections; ///< container for the subsampled output gradient directions GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions BValueMap m_B_ValueMap; double m_B_Value; double m_TargetB_Value; std::vector m_WeightsVector; std::vector > m_ShellInterpolationMatrixVector; std::vector m_bZeroIndicesSplitVectors; IndicesVector m_allDirectionsIndicies; unsigned int m_allDirectionsSize; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkMultiShellAdcAverageReconstructionImageFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h index ea1bc9b0eb..99bbd561a8 100644 --- a/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h +++ b/Modules/DiffusionImaging/DiffusionCore/IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h @@ -1,123 +1,123 @@ /*=================================================================== 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 __mitkDiffusionImage__h #define __mitkDiffusionImage__h #include "mitkImage.h" #include "itkVectorImage.h" #include "itkVectorImageToImageAdaptor.h" #include #include namespace mitk { /** * \brief this class encapsulates diffusion volumes (vectorimages not * yet supported by mitkImage) */ template class DiffusionImage : public Image { public: typedef TPixelType PixelType; typedef typename itk::VectorImage ImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef itk::VectorImageToImageAdaptor< TPixelType, 3 > AdaptorType; typedef vnl_matrix_fixed< double, 3, 3 > MeasurementFrameType; // BValue Map // key := b-Value // value := indicesVector (containing corresponding gradient directions for a b-Value-Shell typedef std::vector< unsigned int > IndicesVector; - typedef std::map< double , IndicesVector > BValueMap; + typedef std::map< unsigned int , IndicesVector > BValueMap; mitkClassMacro( DiffusionImage, Image ) itkNewMacro(Self) void AverageRedundantGradients(double precision); GradientDirectionContainerType::Pointer CalcAveragedDirectionSet(double precision, GradientDirectionContainerType::Pointer directions); void CorrectDKFZBrokenGradientScheme(double precision); typename ImageType::Pointer GetVectorImage() { return m_VectorImage; } void SetVectorImage(typename ImageType::Pointer image ) { this->m_VectorImage = image; } void InitializeFromVectorImage(); void SetDisplayIndexForRendering(int displayIndex); GradientDirectionContainerType::Pointer GetDirectionsWithoutMeasurementFrame() { return m_OriginalDirections; } GradientDirectionContainerType::Pointer GetDirections() { return m_Directions; } void SetDirections( GradientDirectionContainerType::Pointer directions ) { this->m_OriginalDirections = directions; ApplyMeasurementFrame(); } void SetDirections(const std::vector > directions); MeasurementFrameType GetMeasurementFrame() { return m_MeasurementFrame; } void SetMeasurementFrame( MeasurementFrameType mFrame ) { this->m_MeasurementFrame = mFrame; this->ApplyMeasurementFrame(); } bool AreAlike(GradientDirectionType g1, GradientDirectionType g2, double precision); int GetNumDirections(); int GetNumB0(); float GetB_Value(int i); bool IsMultiBval(); void UpdateBValueMap(); IndicesVector GetB0Indices(); itkGetMacro(B_Value, float) itkSetMacro(B_Value, float) BValueMap GetB_ValueMap(){ return m_B_ValueMap; } void AddDirectionsContainerObserver(); void RemoveDirectionsContainerObserver(); protected: DiffusionImage(); virtual ~DiffusionImage(); void ApplyMeasurementFrame(); typename ImageType::Pointer m_VectorImage; GradientDirectionContainerType::Pointer m_Directions; GradientDirectionContainerType::Pointer m_OriginalDirections; float m_B_Value; typename AdaptorType::Pointer m_VectorImageAdaptor; int m_DisplayIndex; MeasurementFrameType m_MeasurementFrame; BValueMap m_B_ValueMap; unsigned long m_DirectionsObserverTag; }; } // namespace mitk #include "mitkDiffusionImage.txx" #endif /* __mitkDiffusionImage__h */ diff --git a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp index 2e0f65c932..e7fe6c5b5c 100644 --- a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp @@ -1,250 +1,250 @@ /*=================================================================== 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 "mitkDiffusionFunctionCollection.h" #include #include "mitkVector.h" // for Windows #ifndef M_PI #define M_PI 3.14159265358979323846 #endif // Namespace ::SH #include #include #include // Namespace ::Gradients #include "itkVectorContainer.h" #include "vnl/vnl_vector.h" //------------------------- SH-function ------------------------------------ double mitk::sh::factorial(int number) { if(number <= 1) return 1; double result = 1.0; for(int i=1; i<=number; i++) result *= i; return result; } void mitk::sh::Cart2Sph(double x, double y, double z, double *cart) { double phi, th, rad; rad = sqrt(x*x+y*y+z*z); if( rad < mitk::eps ) { th = M_PI/2; phi = M_PI/2; } else { th = acos(z/rad); phi = atan2(y, x); } cart[0] = phi; cart[1] = th; cart[2] = rad; } double mitk::sh::legendre0(int l) { if( l%2 != 0 ) { return 0; } else { double prod1 = 1.0; for(int i=1;i mitk::gradients::GetAllUniqueDirections(const std::map > & refBValueMap, const GradientDirectionContainerType * refGradientsContainer ) +std::vector mitk::gradients::GetAllUniqueDirections(const BValueMap & refBValueMap, GradientDirectionContainerType *refGradientsContainer ) { IndiciesVector directioncontainer; BValueMap::const_iterator mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); IndiciesVector::iterator containerIt = directioncontainer.begin(); bool directionExist = false; while(containerIt != directioncontainer.end()) { if (fabs(dot(refGradientsContainer->ElementAt(*containerIt), refGradientsContainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { directioncontainer.push_back(wntIndex); } } } return directioncontainer; } -bool mitk::gradients::CheckForDifferingShellDirections(const std::map > & refBValueMap, const GradientDirectionContainerType * refGradientsContainer) +bool mitk::gradients::CheckForDifferingShellDirections(const BValueMap & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer) { BValueMap::const_iterator mapIterator = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator++; //skip bzero Values for( ; mapIterator != refBValueMap.end(); mapIterator++){ BValueMap::const_iterator mapIterator_2 = refBValueMap.begin(); if(refBValueMap.find(0) != refBValueMap.end() && refBValueMap.size() > 1) mapIterator_2++; //skip bzero Values for( ; mapIterator_2 != refBValueMap.end(); mapIterator_2++){ if(mapIterator_2 == mapIterator) continue; IndiciesVector currentShell = mapIterator->second; IndiciesVector testShell = mapIterator_2->second; for (unsigned int i = 0; i< currentShell.size(); i++) if (fabs(dot(refGradientsContainer->ElementAt(currentShell[i]), refGradientsContainer->ElementAt(testShell[i]))) <= 0.9998) { return true; } } } return false; } template double mitk::gradients::dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ) { double result = (v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2]) / (v1.two_norm() * v2.two_norm()); return result ; } vnl_matrix mitk::gradients::ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer) { vnl_matrix Q(3, refShell.size()); for(unsigned int i = 0; i < refShell.size(); i++) { GradientDirectionType dir = refGradientsContainer->ElementAt(refShell[i]); double x = dir.normalize().get(0); double y = dir.normalize().get(1); double z = dir.normalize().get(2); double cart[3]; mitk::sh::Cart2Sph(x,y,z,cart); Q(0,i) = cart[0]; Q(1,i) = cart[1]; Q(2,i) = cart[2]; } return Q; } vnl_matrix mitk::gradients::ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder) { vnl_matrix SHBasisOutput(QBallReference.cols(), (LOrder+1)*(LOrder+2)*0.5); for(unsigned int i=0; i< SHBasisOutput.rows(); i++) for(int k = 0; k <= LOrder; k += 2) for(int m =- k; m <= k; m++) { int j = ( k * k + k + 2 ) / 2 + m - 1; double phi = QBallReference(0,i); double th = QBallReference(1,i); SHBasisOutput(i,j) = mitk::sh::Yj(m,k,th,phi); } return SHBasisOutput; } mitk::gradients::GradientDirectionContainerType::Pointer mitk::gradients::CreateNormalizedUniqueGradientDirectionContainer(const mitk::gradients::BValueMap & bValueMap, const GradientDirectionContainerType *origninalGradentcontainer) { mitk::gradients::GradientDirectionContainerType::Pointer directioncontainer = mitk::gradients::GradientDirectionContainerType::New(); BValueMap::const_iterator mapIterator = bValueMap.begin(); if(bValueMap.find(0) != bValueMap.end() && bValueMap.size() > 1){ mapIterator++; //skip bzero Values vnl_vector_fixed vec; vec.fill(0.0); directioncontainer->push_back(vec); } for( ; mapIterator != bValueMap.end(); mapIterator++){ IndiciesVector currentShell = mapIterator->second; while(currentShell.size()>0) { unsigned int wntIndex = currentShell.back(); currentShell.pop_back(); mitk::gradients::GradientDirectionContainerType::Iterator containerIt = directioncontainer->Begin(); bool directionExist = false; while(containerIt != directioncontainer->End()) { if (fabs(dot(containerIt.Value(), origninalGradentcontainer->ElementAt(wntIndex))) > 0.9998) { directionExist = true; break; } containerIt++; } if(!directionExist) { GradientDirectionType dir(origninalGradentcontainer->ElementAt(wntIndex)); directioncontainer->push_back(dir.normalize()); } } } return directioncontainer; } diff --git a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h index dcfedfd980..cdff3c6495 100644 --- a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h @@ -1,62 +1,62 @@ /*=================================================================== 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 __mitkDiffusionFunctionCollection_h_ #define __mitkDiffusionFunctionCollection_h_ #include "DiffusionCoreExports.h" #include "vnl/vnl_vector.h" #include "vnl/vnl_vector_fixed.h" #include "itkVectorContainer.h" namespace mitk{ class DiffusionCore_EXPORT sh { public: static double factorial(int number); static void Cart2Sph(double x, double y, double z, double* cart); static double legendre0(int l); static double spherical_harmonic(int m,int l,double theta,double phi, bool complexPart); static double Yj(int m, int k, double theta, double phi); }; class DiffusionCore_EXPORT gradients { private: typedef std::vector IndiciesVector; - typedef std::map BValueMap; - typedef vnl_vector_fixed< double, 3 > GradientDirectionType; - typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; + typedef std::map BValueMap; + typedef itk::VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType; public: - static std::vector GetAllUniqueDirections(const std::map > & refBValueMap, const GradientDirectionContainerType *refGradientsContainer ); - static bool CheckForDifferingShellDirections(const std::map > & refBValueMap, const GradientDirectionContainerType * refGradientsContainer); + static std::vector GetAllUniqueDirections(const BValueMap &bValueMap, GradientDirectionContainerType *refGradientsContainer ); + + static bool CheckForDifferingShellDirections(const BValueMap &bValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer); static vnl_matrix ComputeSphericalHarmonicsBasis(const vnl_matrix & QBallReference, const unsigned int & LOrder); static vnl_matrix ComputeSphericalFromCartesian(const IndiciesVector & refShell, const GradientDirectionContainerType * refGradientsContainer); static mitk::gradients::GradientDirectionContainerType::Pointer CreateNormalizedUniqueGradientDirectionContainer(const BValueMap &bValueMap, const GradientDirectionContainerType * origninalGradentcontainer); template static double dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ); }; } #endif //__mitkDiffusionFunctionCollection_h_ diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h index 8103972b57..828e1cb389 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h @@ -1,122 +1,122 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_h_ #define _itk_ElectrostaticRepulsionDiffusionGradientReductionFilter_h_ #include #include #include namespace itk { /** * \brief Select subset of the input vectors equally distributed over the sphere using an iterative electrostatic repulsion strategy. */ template class ElectrostaticRepulsionDiffusionGradientReductionFilter : public ImageToImageFilter, itk::VectorImage > { public: typedef ElectrostaticRepulsionDiffusionGradientReductionFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< itk::VectorImage, itk::VectorImage > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(ElectrostaticRepulsionDiffusionGradientReductionFilter, ImageToImageFilter) typedef TInputScalarType InputScalarType; typedef itk::VectorImage InputImageType; typedef typename InputImageType::PixelType InputPixelType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef std::vector IndicesVector; - typedef std::map BValueMap; + typedef std::map BValueMap; itkGetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkSetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkGetMacro(GradientDirections, GradientDirectionContainerType::Pointer) itkSetMacro(GradientDirections, GradientDirectionContainerType::Pointer) IndicesVector GetUsedGradientIndices(){return m_UsedGradientIndices;} void SetOriginalBValueMap(BValueMap inp){m_OriginalBValueMap = inp;} void SetShellSelectionBValueMap(BValueMap inp){m_InputBValueMap = inp;} void SetNumGradientDirections(std::vector numDirs){m_NumGradientDirections = numDirs;} protected: ElectrostaticRepulsionDiffusionGradientReductionFilter(); ~ElectrostaticRepulsionDiffusionGradientReductionFilter() {} void GenerateData(); double Costs(); ///< calculates electrostatic energy of current direction set GradientDirectionContainerType::Pointer m_GradientDirections; ///< container for the subsampled output gradient directions GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions IndicesVector m_UsedGradientIndices; IndicesVector m_UnusedGradientIndices; IndicesVector m_BaselineImageIndices; BValueMap m_OriginalBValueMap; BValueMap m_InputBValueMap; std::vector m_NumGradientDirections; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx" #endif #endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp index 1576d6c260..f1e7fba451 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp @@ -1,1062 +1,1062 @@ /*=================================================================== 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. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "QmitkQBallReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include // itk includes #include "itkTimeProbe.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "itkDiffusionQballReconstructionImageFilter.h" #include "itkAnalyticalDiffusionQballReconstructionImageFilter.h" #include "itkDiffusionMultiShellQballReconstructionImageFilter.h" #include "itkVectorContainer.h" #include "mitkQBallImage.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "berryIStructuredSelection.h" #include "berryIWorkbenchWindow.h" #include "berryISelectionService.h" #include const std::string QmitkQBallReconstructionView::VIEW_ID = "org.mitk.views.qballreconstruction"; typedef float TTensorPixelType; const int QmitkQBallReconstructionView::nrconvkernels = 252; struct QbrShellSelection { QmitkQBallReconstructionView* m_View; mitk::DataNode * m_Node; std::string m_NodeName; std::vector m_CheckBoxes; QLabel * m_Label; mitk::DiffusionImage * m_Image; typedef mitk::DiffusionImage::BValueMap BValueMap; QbrShellSelection(QmitkQBallReconstructionView* view, mitk::DataNode * node) : m_View(view), m_Node(node), m_NodeName(node->GetName()) { m_Image = dynamic_cast * > (node->GetData()); if(!m_Image){MITK_INFO << "QmitkQBallReconstructionView::QbrShellSelection : fail to initialize DiffusionImage "; return;} GenerateCheckboxes(); } void GenerateCheckboxes() { BValueMap origMap = m_Image->GetB_ValueMap(); BValueMap::iterator itStart = origMap.begin(); itStart++; BValueMap::iterator itEnd = origMap.end(); m_Label = new QLabel(m_NodeName.c_str()); m_Label->setVisible(true); m_View->m_Controls->m_QBallSelectionBox->layout()->addWidget(m_Label); for(BValueMap::iterator it = itStart ; it!= itEnd; it++) { QCheckBox * box = new QCheckBox(QString::number(it->first)); m_View->m_Controls->m_QBallSelectionBox->layout()->addWidget(box); box->setChecked(true); box->setCheckable(true); // box->setVisible(true); m_CheckBoxes.push_back(box); } } void SetVisible(bool vis) { foreach(QCheckBox * box, m_CheckBoxes) { box->setVisible(vis); } } BValueMap GetBValueSelctionMap() { BValueMap inputMap = m_Image->GetB_ValueMap(); BValueMap outputMap; - double val = 0; + unsigned int val = 0; if(inputMap.find(0) == inputMap.end()){ MITK_INFO << "QbrShellSelection: return empty BValueMap from GUI Selection"; return outputMap; }else{ outputMap[val] = inputMap[val]; MITK_INFO << val; } foreach(QCheckBox * box, m_CheckBoxes) { if(box->isChecked()){ val = box->text().toDouble(); outputMap[val] = inputMap[val]; MITK_INFO << val; } } return outputMap; } ~QbrShellSelection() { m_View->m_Controls->m_QBallSelectionBox->layout()->removeWidget(m_Label); delete m_Label; for(std::vector::iterator it = m_CheckBoxes.begin() ; it!= m_CheckBoxes.end(); it++) { m_View->m_Controls->m_QBallSelectionBox->layout()->removeWidget((*it)); delete (*it); } m_CheckBoxes.clear(); } }; using namespace berry; struct QbrSelListener : ISelectionListener { berryObjectMacro(QbrSelListener); QbrSelListener(QmitkQBallReconstructionView* view) { m_View = view; } void DoSelectionChanged(ISelection::ConstPointer selection) { // save current selection in member variable m_View->m_CurrentSelection = selection.Cast(); // do something with the selected items if(m_View->m_CurrentSelection) { bool foundDwiVolume = false; m_View->m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_View->m_Controls->m_InputData->setTitle("Please Select Input Data"); QString selected_images = ""; mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); int at = 0; // iterate selection for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin(); i != m_View->m_CurrentSelection->End(); ++i) { // extract datatree node if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); mitk::DiffusionImage* diffusionImage; // only look at interesting types if(diffusionImage = dynamic_cast * >(node->GetData())) { foundDwiVolume = true; selected_images += QString(node->GetName().c_str()); if(i + 1 != m_View->m_CurrentSelection->End()) selected_images += "\n"; set->InsertElement(at++, node); } } } m_View->GenerateShellSelectionUI(set); m_View->m_Controls->m_DiffusionImageLabel->setText(selected_images); m_View->m_Controls->m_ButtonStandard->setEnabled(foundDwiVolume); if (foundDwiVolume) m_View->m_Controls->m_InputData->setTitle("Input Data"); else m_View->m_Controls->m_DiffusionImageLabel->setText("mandatory"); } } void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection) { // check, if selection comes from datamanager if (part) { QString partname(part->GetPartName().c_str()); if(partname.compare("Datamanager")==0) { // apply selection DoSelectionChanged(selection); } } } QmitkQBallReconstructionView* m_View; }; // --------------- QmitkQBallReconstructionView----------------- // QmitkQBallReconstructionView::QmitkQBallReconstructionView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL) { } QmitkQBallReconstructionView::QmitkQBallReconstructionView(const QmitkQBallReconstructionView& other) { Q_UNUSED(other); throw std::runtime_error("Copy constructor not implemented"); } QmitkQBallReconstructionView::~QmitkQBallReconstructionView() { this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->RemovePostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener); } void QmitkQBallReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkQBallReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); m_Controls->m_DiffusionImageLabel->setText("mandatory"); QStringList items; items << "2" << "4" << "6" << "8" << "10" << "12"; m_Controls->m_QBallReconstructionMaxLLevelComboBox->addItems(items); m_Controls->m_QBallReconstructionMaxLLevelComboBox->setCurrentIndex(1); MethodChoosen(m_Controls->m_QBallReconstructionMethodComboBox->currentIndex()); #ifndef DIFFUSION_IMAGING_EXTENDED m_Controls->m_QBallReconstructionMethodComboBox->removeItem(3); #endif AdvancedCheckboxClicked(); } m_SelListener = berry::ISelectionListener::Pointer(new QbrSelListener(this)); this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->AddPostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener); berry::ISelection::ConstPointer sel( this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); m_CurrentSelection = sel.Cast(); m_SelListener.Cast()->DoSelectionChanged(sel); } void QmitkQBallReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkQBallReconstructionView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkQBallReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_ButtonStandard), SIGNAL(clicked()), this, SLOT(ReconstructStandard()) ); connect( (QObject*)(m_Controls->m_AdvancedCheckbox), SIGNAL(clicked()), this, SLOT(AdvancedCheckboxClicked()) ); connect( (QObject*)(m_Controls->m_QBallReconstructionMethodComboBox), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) ); } } void QmitkQBallReconstructionView::OnSelectionChanged( std::vector nodes ) { } void QmitkQBallReconstructionView::Activated() { QmitkFunctionality::Activated(); berry::ISelection::ConstPointer sel( this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); m_CurrentSelection = sel.Cast(); m_SelListener.Cast()->DoSelectionChanged(sel); } void QmitkQBallReconstructionView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkQBallReconstructionView::ReconstructStandard() { int index = m_Controls->m_QBallReconstructionMethodComboBox->currentIndex(); #ifndef DIFFUSION_IMAGING_EXTENDED if(index>=3) { index = index + 1; } #endif switch(index) { case 0: { // Numerical Reconstruct(0,0); break; } case 1: { // Standard Reconstruct(1,0); break; } case 2: { // Solid Angle Reconstruct(1,6); break; } case 3: { // Constrained Solid Angle Reconstruct(1,7); break; } case 4: { // ADC Reconstruct(1,4); break; } case 5: { // Raw Signal Reconstruct(1,5); break; } case 6: { // Q-Ball reconstruction Reconstruct(2,0); break; } } } void QmitkQBallReconstructionView::MethodChoosen(int method) { #ifndef DIFFUSION_IMAGING_EXTENDED if(method>=3) { method = method + 1; } #endif m_Controls->m_QBallSelectionBox->setHidden(true); m_Controls->m_OutputCoeffsImage->setHidden(true); if (method==0) m_Controls->m_ShFrame->setVisible(false); else m_Controls->m_ShFrame->setVisible(true); switch(method) { case 0: m_Controls->m_Description->setText("Numerical recon. (Tuch 2004)"); break; case 1: m_Controls->m_Description->setText("Spherical harmonics recon. (Descoteaux 2007)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 2: m_Controls->m_Description->setText("SH recon. with solid angle consideration (Aganj 2009)"); m_Controls->m_OutputCoeffsImage->setHidden(false); break; case 3: m_Controls->m_Description->setText("SH solid angle with non-neg. constraint (Goh 2009)"); break; case 4: m_Controls->m_Description->setText("SH recon. of the plain ADC-profiles"); break; case 5: m_Controls->m_Description->setText("SH recon. of the raw diffusion signal"); break; case 6: m_Controls->m_Description->setText("SH recon. of the multi shell diffusion signal (Aganj 2010)"); m_Controls->m_QBallSelectionBox->setHidden(false); m_Controls->m_OutputCoeffsImage->setHidden(false); break; } } void QmitkQBallReconstructionView::AdvancedCheckboxClicked() { bool check = m_Controls->m_AdvancedCheckbox->isChecked(); m_Controls->m_QBallReconstructionMaxLLevelTextLabel_2->setVisible(check); m_Controls->m_QBallReconstructionMaxLLevelComboBox->setVisible(check); m_Controls->m_QBallReconstructionLambdaTextLabel_2->setVisible(check); m_Controls->m_QBallReconstructionLambdaLineEdit->setVisible(check); m_Controls->m_QBallReconstructionThresholdLabel_2->setVisible(check); m_Controls->m_QBallReconstructionThreasholdEdit->setVisible(check); m_Controls->label_2->setVisible(check); m_Controls->frame_2->setVisible(check); } void QmitkQBallReconstructionView::Reconstruct(int method, int normalization) { if (m_CurrentSelection) { mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); int at = 0; for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); i != m_CurrentSelection->End(); ++i) { if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { set->InsertElement(at++, node); } } } if(method == 0) { NumericalQBallReconstruction(set, normalization); } else { #if BOOST_VERSION / 100000 > 0 #if BOOST_VERSION / 100 % 1000 > 34 if(method == 1) { AnalyticalQBallReconstruction(set, normalization); } if(method == 2) { MultiQBallReconstruction(set); } #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; QMessageBox::warning(NULL,"ERROR","Boost 1.35 minimum required"); #endif #else std::cout << "ERROR: Boost 1.35 minimum required" << std::endl; QMessageBox::warning(NULL,"ERROR","Boost 1.35 minimum required"); #endif } } } void QmitkQBallReconstructionView::NumericalQBallReconstruction (mitk::DataStorage::SetOfObjects::Pointer inImages, int normalization) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // QBALL RECONSTRUCTION clock.Start(); MITK_INFO << "QBall reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "QBall reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::DiffusionQballReconstructionImageFilter QballReconstructionImageFilterType; QballReconstructionImageFilterType::Pointer filter = QballReconstructionImageFilterType::New(); filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); filter->SetBValue(vols->GetB_Value()); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); switch(normalization) { case 0: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); break; } case 1: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO_B_VALUE); break; } case 2: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO); break; } case 3: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_NONE); break; } default: { filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD); } } filter->Update(); clock.Stop(); MITK_DEBUG << "took " << clock.GetMeanTime() << "s." ; // ODFs TO DATATREE mitk::QBallImage::Pointer image = mitk::QBallImage::New(); image->InitializeByItk( filter->GetOutput() ); //image->SetImportVolume( filter->GetOutput()->GetBufferPointer(), 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_QN%1").arg(normalization); SetDefaultNodeProperties(node, newname.toStdString()); nodes.push_back(node); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return ; } } void QmitkQBallReconstructionView::AnalyticalQBallReconstruction( mitk::DataStorage::SetOfObjects::Pointer inImages, int normalization) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; std::vector lambdas; float minLambda = m_Controls->m_QBallReconstructionLambdaLineEdit->value(); lambdas.push_back(minLambda); int nLambdas = lambdas.size(); QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles*nLambdas); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector* nodes = new std::vector(); while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name",nodename); itemiter++; // QBALL RECONSTRUCTION clock.Start(); MITK_INFO << "QBall reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "QBall reconstruction for %s", nodename.c_str()).toAscii()); for(int i=0; im_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedAnalyticalQBallReconstruction<2>(vols, currentLambda, nodename, nodes, normalization); break; } case 1: { TemplatedAnalyticalQBallReconstruction<4>(vols, currentLambda, nodename, nodes, normalization); break; } case 2: { TemplatedAnalyticalQBallReconstruction<6>(vols, currentLambda, nodename, nodes, normalization); break; } case 3: { TemplatedAnalyticalQBallReconstruction<8>(vols, currentLambda, nodename, nodes, normalization); break; } case 4: { TemplatedAnalyticalQBallReconstruction<10>(vols, currentLambda, nodename, nodes, normalization); break; } case 5: { TemplatedAnalyticalQBallReconstruction<12>(vols, currentLambda, nodename, nodes, normalization); break; } } clock.Stop(); MITK_DEBUG << "took " << clock.GetMeanTime() << "s." ; mitk::ProgressBar::GetInstance()->Progress(); } } std::vector::iterator nodeIt; for(nodeIt = nodes->begin(); nodeIt != nodes->end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); m_MultiWidget->RequestUpdate(); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return; } } template void QmitkQBallReconstructionView::TemplatedAnalyticalQBallReconstruction( mitk::DiffusionImage* vols, float lambda, std::string nodename, std::vector* nodes, int normalization) { typedef itk::AnalyticalDiffusionQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); filter->SetBValue(vols->GetB_Value()); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); filter->SetLambda(lambda); switch(normalization) { case 0: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); break; } case 1: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO_B_VALUE); break; } case 2: { filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO); break; } case 3: { filter->SetNormalizationMethod(FilterType::QBAR_NONE); break; } case 4: { filter->SetNormalizationMethod(FilterType::QBAR_ADC_ONLY); break; } case 5: { filter->SetNormalizationMethod(FilterType::QBAR_RAW_SIGNAL); break; } case 6: { filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE); break; } case 7: { filter->SetNormalizationMethod(FilterType::QBAR_NONNEG_SOLID_ANGLE); break; } default: { filter->SetNormalizationMethod(FilterType::QBAR_STANDARD); } } filter->Update(); // ODFs TO DATATREE mitk::QBallImage::Pointer image = mitk::QBallImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_QA%1").arg(normalization); SetDefaultNodeProperties(node, newname.toStdString()); nodes->push_back(node); if(m_Controls->m_OutputCoeffsImage->isChecked()) { mitk::Image::Pointer coeffsImage = mitk::Image::New(); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); mitk::DataNode::Pointer coeffsNode=mitk::DataNode::New(); coeffsNode->SetData( coeffsImage ); coeffsNode->SetProperty( "name", mitk::StringProperty::New( QString(nodename.c_str()).append("_coeffs").toStdString()) ); coeffsNode->SetVisibility(false); nodes->push_back(coeffsNode); } } void QmitkQBallReconstructionView::MultiQBallReconstruction( mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; std::vector lambdas; float minLambda = m_Controls->m_QBallReconstructionLambdaLineEdit->value(); lambdas.push_back(minLambda); int nLambdas = lambdas.size(); QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles*nLambdas); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector* nodes = new std::vector(); while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); const mitk::DataNode * nodePointer = (*itemiter).GetPointer(); std::string nodename; (*itemiter)->GetStringProperty("name",nodename); itemiter++; // QBALL RECONSTRUCTION clock.Start(); MITK_INFO << "QBall reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "QBall reconstruction for %s", nodename.c_str()).toAscii()); for(int i=0; im_QBallReconstructionMaxLLevelComboBox->currentIndex()) { case 0: { TemplatedMultiQBallReconstruction<2>(vols, currentLambda, nodePointer, nodes); break; } case 1: { TemplatedMultiQBallReconstruction<4>(vols, currentLambda, nodePointer, nodes); break; } case 2: { TemplatedMultiQBallReconstruction<6>(vols, currentLambda, nodePointer, nodes); break; } case 3: { TemplatedMultiQBallReconstruction<8>(vols, currentLambda, nodePointer, nodes); break; } case 4: { TemplatedMultiQBallReconstruction<10>(vols, currentLambda, nodePointer, nodes); break; } case 5: { TemplatedMultiQBallReconstruction<12>(vols, currentLambda, nodePointer, nodes); break; } } clock.Stop(); MITK_DEBUG << "took " << clock.GetMeanTime() << "s." ; mitk::ProgressBar::GetInstance()->Progress(); } } std::vector::iterator nodeIt; for(nodeIt = nodes->begin(); nodeIt != nodes->end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); m_MultiWidget->RequestUpdate(); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); } catch (itk::ExceptionObject &ex) { MITK_INFO << ex ; QMessageBox::information(0, "Reconstruction not possible:", ex.GetDescription()); return ; } } template void QmitkQBallReconstructionView::TemplatedMultiQBallReconstruction( mitk::DiffusionImage* vols, float lambda, const mitk::DataNode * dataNodePointer, std::vector* nodes) { typedef itk::DiffusionMultiShellQballReconstructionImageFilter FilterType; typename FilterType::Pointer filter = FilterType::New(); std::string nodename; dataNodePointer->GetStringProperty("name",nodename); filter->SetBValueMap(m_ShellSelectorMap[dataNodePointer]->GetBValueSelctionMap()); filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage(), vols->GetB_Value() ); filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->value() ); filter->SetLambda(lambda); filter->Update(); // ODFs TO DATATREE mitk::QBallImage::Pointer image = mitk::QBallImage::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_QAMultiShell"); SetDefaultNodeProperties(node, newname.toStdString()); nodes->push_back(node); if(m_Controls->m_OutputCoeffsImage->isChecked()) { mitk::Image::Pointer coeffsImage = mitk::Image::New(); coeffsImage->InitializeByItk( filter->GetCoefficientImage().GetPointer() ); coeffsImage->SetVolume( filter->GetCoefficientImage()->GetBufferPointer() ); mitk::DataNode::Pointer coeffsNode=mitk::DataNode::New(); coeffsNode->SetData( coeffsImage ); coeffsNode->SetProperty( "name", mitk::StringProperty::New( QString(nodename.c_str()).append("_coeffs").toStdString()) ); nodes->push_back(coeffsNode); } } void QmitkQBallReconstructionView::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( "opacity", mitk::FloatProperty::New(1.0f) ); node->SetProperty( "name", mitk::StringProperty::New(name) ); } //node->SetProperty( "volumerendering", mitk::BoolProperty::New( false ) ); //node->SetProperty( "use color", mitk::BoolProperty::New( true ) ); //node->SetProperty( "texture interpolation", mitk::BoolProperty::New( true ) ); //node->SetProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() ); //node->SetProperty( "layer", mitk::IntProperty::New(0)); //node->SetProperty( "in plane resample extent by geometry", mitk::BoolProperty::New( false ) ); //node->SetOpacity(1.0f); //node->SetColor(1.0,1.0,1.0); //node->SetVisibility(true); //node->SetProperty( "IsQBallVolume", mitk::BoolProperty::New( true ) ); //mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); //mitk::LevelWindow levelwindow; //// levelwindow.SetAuto( image ); //levWinProp->SetLevelWindow( levelwindow ); //node->GetPropertyList()->SetPropertx( "levelwindow", levWinProp ); //// add a default rainbow lookup table for color mapping //if(!node->GetProperty("LookupTable")) //{ // mitk::LookupTable::Pointer mitkLut = mitk::LookupTable::New(); // vtkLookupTable* vtkLut = mitkLut->GetVtkLookupTable(); // vtkLut->SetHueRange(0.6667, 0.0); // vtkLut->SetTableRange(0.0, 20.0); // vtkLut->Build(); // mitk::LookupTableProperty::Pointer mitkLutProp = mitk::LookupTableProperty::New(); // mitkLutProp->SetLookupTable(mitkLut); // node->SetProperty( "LookupTable", mitkLutProp ); //} //if(!node->GetProperty("binary")) // node->SetProperty( "binary", mitk::BoolProperty::New( false ) ); //// add a default transfer function //mitk::TransferFunction::Pointer tf = mitk::TransferFunction::New(); //node->SetProperty ( "TransferFunction", mitk::TransferFunctionProperty::New ( tf.GetPointer() ) ); //// set foldername as string property //mitk::StringProperty::Pointer nameProp = mitk::StringProperty::New( name ); //node->SetProperty( "name", nameProp ); void QmitkQBallReconstructionView::GenerateShellSelectionUI(mitk::DataStorage::SetOfObjects::Pointer set) { std::map tempMap; const mitk::DataStorage::SetOfObjects::iterator setEnd( set->end() ); mitk::DataStorage::SetOfObjects::iterator NodeIt( set->begin() ); while(NodeIt != setEnd) { if(m_ShellSelectorMap.find( (*NodeIt).GetPointer() ) != m_ShellSelectorMap.end()) { tempMap[(*NodeIt).GetPointer()] = m_ShellSelectorMap[(*NodeIt).GetPointer()]; m_ShellSelectorMap.erase((*NodeIt).GetPointer()); }else { tempMap[(*NodeIt).GetPointer()] = new QbrShellSelection(this, (*NodeIt) ); tempMap[(*NodeIt).GetPointer()]->SetVisible(true); } NodeIt++; } for(std::map::iterator it = m_ShellSelectorMap.begin(); it != m_ShellSelectorMap.end();it ++) { delete it->second; } m_ShellSelectorMap.clear(); m_ShellSelectorMap = tempMap; }