diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp index 307019dd83..f3e5d43528 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp @@ -1,150 +1,194 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.txx $ Language: C++ Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ Version: $Revision: 68 $ Copyright (c) INRIA 2010. All rights reserved. See LICENSE.txt for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef _itk_MultiShellAdcAverageReconstructionImageFilter_cpp_ #define _itk_MultiShellAdcAverageReconstructionImageFilter_cpp_ #endif #define _USE_MATH_DEFINES #include "itkMultiShellAdcAverageReconstructionImageFilter.h" #include #include #include #include +#include "mitkDiffusionFunctionCollection.h" namespace itk { template MultiShellAdcAverageReconstructionImageFilter -::MultiShellAdcAverageReconstructionImageFilter() +::MultiShellAdcAverageReconstructionImageFilter(): + m_Interpolation(false) { this->SetNumberOfRequiredInputs( 1 ); } template void MultiShellAdcAverageReconstructionImageFilter ::BeforeThreadedGenerateData() { - // test whether interpolation is necessary - // - Gradeint directions on different shells are different - //# if interpolation is neccesary and any #ShellDirection < 15 --> itkException (No interpolation possible) // test whether BvalueMap contains all necessary information - //# BValueMapSize == 0 --> itkWarning (take BValueMap from inputImage) + if(m_BValueMap.size() == 0) + { + itkWarningMacro(<< "No BValueMap given: create one using GradientDirectionContainer"); + + GradientDirectionContainerType::ConstIterator gdcit; + for( gdcit = m_OriginalGradientDirections->Begin(); gdcit != m_OriginalGradientDirections->End(); ++gdcit) + { + double bValueKey = int(((m_BValue * gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; + m_BValueMap[bValueKey].push_back(gdcit.Index()); + } + } + //# BValueMap contains no bZero --> itkException + if(m_BValueMap.find(0.0) == m_BValueMap.end()) + { + MITK_INFO << "No ReferenceSignal (BZeroImages) found!"; + itkExceptionMacro(<< "No ReferenceSignal (BZeroImages) found!"); + } + + // test whether interpolation is necessary + // - Gradeint directions on different shells are different + m_Interpolation = mitk::gradients::CheckForDifferingShellDirections(m_BValueMap, m_OriginalGradientDirections.GetPointer()); // if INTERPOLATION necessary - // [allDirectionsContainer] Gradient DirectionContainer containing all unique directions - // [sizeAllDirections] size of GradientContainer cointaining all unique directions - /* for each shell - * - calculate maxShOrder - * - calculate Weights [Weigthing = shell_size / max_shell_size] - * - get TragetSHBasis using allDirectionsContainer - * - get ShellSHBasis using currentShellDirections - * - calculate interpolationSHBasis [TargetSHBasis * ShellSHBasis^-1] - * - save interpolationSHBasis - */ + if(m_Interpolation) + { + for(BValueMap::const_iterator it = m_BValueMap.begin();it != m_BValueMap.end(); it++) + { + if((*it).first == 0.0) continue; + // if any #ShellDirection < 15 --> itkException (No interpolation possible) + if((*it).second.size() < 15){ + MITK_INFO << "Abort: No interpolation possible Shell-" << (*it).first << " has less than 15 directions."; + itkExceptionMacro(<<"No interpolation possible"); + } + } + + // [allDirectionsContainer] Gradient DirectionContainer containing all unique directions + IndicesVector allDirectionsIndicies = mitk::gradients::GetAllUniqueDirections(m_BValueMap, m_OriginalGradientDirections); + + // [sizeAllDirections] size of GradientContainer cointaining all unique directions + const int allDirectionsSize = allDirections.size(); + + /* for each shell + * - calculate maxShOrder + * - calculate Weights [Weigthing = shell_size / max_shell_size] + * - get TragetSHBasis using allDirectionsContainer + * - get ShellSHBasis using currentShellDirections + * - calculate interpolationSHBasis [TargetSHBasis * ShellSHBasis^-1] + * - save interpolationSHBasis + */ + } // calculate average b-Value for target b-Value [bVal_t] // calculate target bZero-Value [b0_t] + MITK_INFO << "Input:" << std::endl + << "GradientDirections: " << m_OriginalGradientDirections->Size() << std::endl + << "Shells: " << (m_BValueMap.size() - 1) << std::endl + << "ReferenceImages: " << m_BValueMap.at(0.0).size() << std::endl + << "Interpolation: " << m_Interpolation; + + } template void MultiShellAdcAverageReconstructionImageFilter ::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, int /*threadId*/) { - // Get input gradient image pointer + /* // Get input gradient image pointer typename InputImageType::Pointer inputImage = static_cast< InputImageType * >(ProcessObject::GetInput(0)); // ImageRegionIterator for the input image ImageRegionIterator< InputImageType > iit(inputImage, outputRegionForThread); iit.GoToBegin(); // Get output gradient image pointer typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(ProcessObject::GetOutput(0)); // ImageRegionIterator for the output image ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); const int numShells = m_BValueMap.size()-1; BValueMap::iterator it = m_BValueMap.begin(); //std::vector adcVec = new Vector(numShells); - +*/ // create empty nxm SignalMatrix containing n->signals/directions (in case of interpolation ~ sizeAllDirections otherwise the size of any shell) for m->shells // create nx1 targetSignalVector // ** walking over each Voxel /* for each shell in this Voxel * - get the RawSignal * - interpolate the Signal if necessary using corresponding interpolationSHBasis * - save the (interpolated) ShellSignalVector as the ith column in the SignalMatrix */ // normalize the signals in SignalMatrix [bZeroAverage????] [S/S0 = E] /* for each row in the SignalMatrix * - calculate for each rowentry the ADC [ln(E)/-b = ADC] * - weight the ith ADC entry with the corresponding shellWeighting * - average the Signal over the row [ADC_t] * - calculate target Signal using target b-Value [S_t = b0_t * e^(-bVal_t * ADC_t) * - Save S_t in targetSignalVector */ // outImageIterator set S_t // ** - + /* int vecLength; // initialize output image typename OutputImageType::Pointer outImage = OutputImageType::New(); outImage->SetSpacing( this->GetInput()->GetSpacing() ); outImage->SetOrigin( this->GetInput()->GetOrigin() ); outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction using bZeroDirection+AllDirectionsContainer outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion()); outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() ); outImage->SetVectorLength( vecLength ); // size of 1(bzeroValue) + AllDirectionsContainer outImage->Allocate(); this->SetNumberOfRequiredOutputs (1); this->SetNthOutput (0, outImage); MITK_INFO << "...done"; + */ } } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h index 1eea82d3f3..249a298e63 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h @@ -1,120 +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. ===================================================================*/ -/*========================================================================= -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_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; itkGetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkSetMacro(OriginalGradientDirections, GradientDirectionContainerType::Pointer) itkGetMacro(TargetGradientDirections, GradientDirectionContainerType::Pointer) itkGetMacro(TargetBValue, float) + itkGetMacro(BValue,float) + itkSetMacro(BValue,float) inline void SetOriginalBValueMap(BValueMap inp){m_BValueMap = inp;} protected: MultiShellAdcAverageReconstructionImageFilter(); ~MultiShellAdcAverageReconstructionImageFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int NumberOfThreads ); GradientDirectionContainerType::Pointer m_TargetGradientDirections; ///< 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_BValueMap; + float m_BValue; float m_TargetBValue; + bool m_Interpolation; + }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkMultiShellAdcAverageReconstructionImageFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp index 04c092a739..20902c1114 100644 --- a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.cpp @@ -1,169 +1,169 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "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( std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer ) +std::vector mitk::gradients::GetAllUniqueDirections(const std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer ) { IndiciesVector directioncontainer; - BValueMap::iterator mapIterator = refBValueMap.begin(); + 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(std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer) +bool mitk::gradients::CheckForDifferingShellDirections(const std::map > & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer) { - BValueMap::iterator mapIterator = refBValueMap.begin(); + 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::iterator mapIterator_2 = refBValueMap.begin(); + 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 ; } diff --git a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h index 84103a505e..66f91ffdd1 100644 --- a/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h +++ b/Modules/DiffusionImaging/DiffusionCore/mitkDiffusionFunctionCollection.h @@ -1,61 +1,61 @@ /*=================================================================== 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 "vector.h" #include "map.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); + 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 itk::VectorContainer< unsigned int, vnl_vector_fixed< double, 3 > > GradientDirectionContainerType; public: - static std::vector GetAllUniqueDirections( std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer ); + static std::vector GetAllUniqueDirections(const std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer ); - static bool CheckForDifferingShellDirections(std::map > & refBValueMap, GradientDirectionContainerType *refGradientsContainer); + static bool CheckForDifferingShellDirections(const std::map > & refBValueMap, GradientDirectionContainerType::ConstPointer refGradientsContainer); template static double dot (vnl_vector_fixed< type ,3> const& v1, vnl_vector_fixed< type ,3 > const& v2 ); }; } #endif //__mitkDiffusionFunctionCollection_h_ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp index 85afebbb44..74b210b5e8 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp @@ -1,634 +1,635 @@ /*=================================================================== 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 "QmitkPreprocessingView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include // itk includes #include "itkTimeProbe.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "itkBrainMaskExtractionImageFilter.h" #include "itkCastImageFilter.h" #include "itkVectorContainer.h" #include #include #include // mitk includes #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.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 #include #include const std::string QmitkPreprocessingView::VIEW_ID = "org.mitk.views.preprocessing"; #define DI_INFO MITK_INFO("DiffusionImaging") typedef float TTensorPixelType; QmitkPreprocessingView::QmitkPreprocessingView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL), m_DiffusionImage(NULL) { } QmitkPreprocessingView::QmitkPreprocessingView(const QmitkPreprocessingView& other) { Q_UNUSED(other) throw std::runtime_error("Copy constructor not implemented"); } QmitkPreprocessingView::~QmitkPreprocessingView() { } void QmitkPreprocessingView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkPreprocessingViewControls; m_Controls->setupUi(parent); this->CreateConnections(); m_Controls->m_MeasurementFrameTable->horizontalHeader()->setResizeMode(QHeaderView::Stretch); m_Controls->m_MeasurementFrameTable->verticalHeader()->setResizeMode(QHeaderView::Stretch); } } void QmitkPreprocessingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkPreprocessingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkPreprocessingView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_ButtonAverageGradients), SIGNAL(clicked()), this, SLOT(AverageGradients()) ); connect( (QObject*)(m_Controls->m_ButtonExtractB0), SIGNAL(clicked()), this, SLOT(ExtractB0()) ); connect( (QObject*)(m_Controls->m_ModifyMeasurementFrame), SIGNAL(clicked()), this, SLOT(DoApplyMesurementFrame()) ); connect( (QObject*)(m_Controls->m_ReduceGradientsButton), SIGNAL(clicked()), this, SLOT(DoReduceGradientDirections()) ); connect( (QObject*)(m_Controls->m_ShowGradientsButton), SIGNAL(clicked()), this, SLOT(DoShowGradientDirections()) ); connect( (QObject*)(m_Controls->m_MirrorGradientToHalfSphereButton), SIGNAL(clicked()), this, SLOT(DoHalfSphereGradientDirections()) ); connect( (QObject*)(m_Controls->m_MergeDwisButton), SIGNAL(clicked()), this, SLOT(MergeDwis()) ); connect( (QObject*)(m_Controls->m_AdcAverage), SIGNAL(clicked()), this, SLOT(DoAdcAverage()) ); } } void QmitkPreprocessingView::DoAdcAverage() { if (m_DiffusionImage.IsNull()) return; typedef mitk::DiffusionImage DiffusionImageType; typedef itk::MultiShellAdcAverageReconstructionImageFilter FilterType; typedef DiffusionImageType::BValueMap BValueMap; BValueMap originalShellMap = m_DiffusionImage->GetB_ValueMap(); GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetDirections(); FilterType::Pointer filter = FilterType::New(); filter->SetInput(m_DiffusionImage->GetVectorImage()); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); + filter->SetBValue(m_DiffusionImage->GetB_Value()); filter->Update(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value( filter->GetTargetBValue() ); image->SetDirections( filter->GetTargetGradientDirections() ); image->SetMeasurementFrame( m_DiffusionImage->GetMeasurementFrame() ); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); imageNode->SetName((name+"_averaged").toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); } void QmitkPreprocessingView::OnSelectionChanged( std::vector nodes ) { bool foundDwiVolume = false; m_DiffusionImage = NULL; m_SelectedDiffusionNodes.clear(); // iterate selection for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast*>(node->GetData()) ) { foundDwiVolume = true; m_DiffusionImage = dynamic_cast*>(node->GetData()); m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_SelectedDiffusionNodes.push_back(node); } } m_Controls->m_ButtonAverageGradients->setEnabled(foundDwiVolume); m_Controls->m_ButtonExtractB0->setEnabled(foundDwiVolume); m_Controls->m_CheckExtractAll->setEnabled(foundDwiVolume); m_Controls->m_ModifyMeasurementFrame->setEnabled(foundDwiVolume); m_Controls->m_MeasurementFrameTable->setEnabled(foundDwiVolume); m_Controls->m_ReduceGradientsButton->setEnabled(foundDwiVolume); m_Controls->m_ShowGradientsButton->setEnabled(foundDwiVolume); m_Controls->m_MirrorGradientToHalfSphereButton->setEnabled(foundDwiVolume); m_Controls->m_MergeDwisButton->setEnabled(foundDwiVolume); foreach(QCheckBox * box, m_ReduceGradientCheckboxes) { m_Controls->m_ReductionFrame->layout()->removeWidget(box); delete box; } foreach(QSpinBox * box, m_ReduceGradientSpinboxes) { m_Controls->m_ReductionFrame->layout()->removeWidget(box); delete box; } m_ReduceGradientCheckboxes.clear(); m_ReduceGradientSpinboxes.clear(); if (foundDwiVolume) { m_Controls->m_InputData->setTitle("Input Data"); vnl_matrix_fixed< double, 3, 3 > mf = m_DiffusionImage->GetMeasurementFrame(); for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); delete item; item = new QTableWidgetItem(); item->setTextAlignment(Qt::AlignCenter | Qt::AlignVCenter); item->setText(QString::number(mf.get(r,c))); m_Controls->m_MeasurementFrameTable->setItem(r,c,item); } typedef mitk::DiffusionImage::BValueMap BValueMap; typedef mitk::DiffusionImage::BValueMap::iterator BValueMapIterator; BValueMap bValMap = m_DiffusionImage->GetB_ValueMap(); BValueMapIterator it = bValMap.begin(); m_Controls->m_BvalueTable->clear(); m_Controls->m_BvalueTable->setRowCount(bValMap.size() ); QStringList headerList; headerList << "b-Value" << "Number of gradients"; m_Controls->m_BvalueTable->setHorizontalHeaderLabels(headerList); QCheckBox* checkBox; QSpinBox* spinBox; int i = 0 ; for(;it != bValMap.end(); it++) { m_Controls->m_BvalueTable->setItem(i,0,new QTableWidgetItem(QString::number(it->first))); m_Controls->m_BvalueTable->setItem(i,1,new QTableWidgetItem(QString::number(it->second.size()))); // Reduce Gradients GUI adaption if(it->first != 0 && bValMap.size() > 1){ checkBox = new QCheckBox(QString::number(it->first) + " with " + QString::number(it->second.size()) + " directions"); checkBox->setEnabled(true); checkBox->setChecked(true); checkBox->setCheckable(true); m_ReduceGradientCheckboxes.push_back(checkBox); m_Controls->m_ReductionFrame->layout()->addWidget(checkBox); spinBox = new QSpinBox(); spinBox->setValue(std::ceil((float)it->second.size()/2)); spinBox->setMaximum(it->second.size()); spinBox->setMinimum(0); m_ReduceGradientSpinboxes.push_back(spinBox); m_Controls->m_ReductionFrame->layout()->addWidget(spinBox); } i++; } } else { for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); delete item; item = new QTableWidgetItem(); m_Controls->m_MeasurementFrameTable->setItem(r,c,item); } m_Controls->m_BvalueTable->clear(); m_Controls->m_BvalueTable->setRowCount(1); QStringList headerList; headerList << "b-Value" << "Number of gradients"; m_Controls->m_BvalueTable->setHorizontalHeaderLabels(headerList); m_Controls->m_BvalueTable->setItem(0,0,new QTableWidgetItem("-")); m_Controls->m_BvalueTable->setItem(0,1,new QTableWidgetItem("-")); m_Controls->m_DiffusionImageLabel->setText("mandatory"); m_Controls->m_InputData->setTitle("Please Select Input Data"); } } void QmitkPreprocessingView::Activated() { QmitkFunctionality::Activated(); } void QmitkPreprocessingView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkPreprocessingView::DoHalfSphereGradientDirections() { GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetDirections(); for (int j=0; jSize(); j++) if (gradientContainer->at(j)[0]<0) gradientContainer->at(j) = -gradientContainer->at(j); m_DiffusionImage->SetDirections(gradientContainer); } void QmitkPreprocessingView::DoApplyMesurementFrame() { if (m_DiffusionImage.IsNull()) return; vnl_matrix_fixed< double, 3, 3 > mf; for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); if (!item) return; mf[r][c] = item->text().toDouble(); } m_DiffusionImage->SetMeasurementFrame(mf); } void QmitkPreprocessingView::DoShowGradientDirections() { if (m_DiffusionImage.IsNull()) return; int maxIndex = 0; int maxSize = m_DiffusionImage->GetDimension(0); if (maxSizeGetDimension(1)) { maxSize = m_DiffusionImage->GetDimension(1); maxIndex = 1; } if (maxSizeGetDimension(2)) { maxSize = m_DiffusionImage->GetDimension(2); maxIndex = 2; } mitk::Point3D origin = m_DiffusionImage->GetGeometry()->GetOrigin(); mitk::PointSet::Pointer originSet = mitk::PointSet::New(); typedef mitk::DiffusionImage::BValueMap BValueMap; typedef mitk::DiffusionImage::BValueMap::iterator BValueMapIterator; BValueMap bValMap = m_DiffusionImage->GetB_ValueMap(); GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetDirections(); mitk::Geometry3D::Pointer geometry = m_DiffusionImage->GetGeometry(); int shellCount = 1; for(BValueMapIterator it = bValMap.begin(); it!=bValMap.end(); ++it) { mitk::PointSet::Pointer pointset = mitk::PointSet::New(); for (int j=0; jsecond.size(); j++) { mitk::Point3D ip; vnl_vector_fixed< double, 3 > v = gradientContainer->at(it->second[j]); if (v.magnitude()>mitk::eps) { ip[0] = v[0]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[0]-0.5*geometry->GetSpacing()[0] + geometry->GetSpacing()[0]*m_DiffusionImage->GetDimension(0)/2; ip[1] = v[1]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[1]-0.5*geometry->GetSpacing()[1] + geometry->GetSpacing()[1]*m_DiffusionImage->GetDimension(1)/2; ip[2] = v[2]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[2]-0.5*geometry->GetSpacing()[2] + geometry->GetSpacing()[2]*m_DiffusionImage->GetDimension(2)/2; pointset->InsertPoint(j, ip); } else if (originSet->IsEmpty()) { ip[0] = v[0]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[0]-0.5*geometry->GetSpacing()[0] + geometry->GetSpacing()[0]*m_DiffusionImage->GetDimension(0)/2; ip[1] = v[1]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[1]-0.5*geometry->GetSpacing()[1] + geometry->GetSpacing()[1]*m_DiffusionImage->GetDimension(1)/2; ip[2] = v[2]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[2]-0.5*geometry->GetSpacing()[2] + geometry->GetSpacing()[2]*m_DiffusionImage->GetDimension(2)/2; originSet->InsertPoint(j, ip); } } if (it->firstSetData(pointset); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); name += "_Shell_"; name += QString::number(it->first); node->SetName(name.toStdString().c_str()); node->SetProperty("pointsize", mitk::FloatProperty::New((float)maxSize/50)); int b0 = shellCount%2; int b1 = 0; int b2 = 0; if (shellCount>4) b2 = 1; if (shellCount%4 >= 2) b1 = 1; node->SetProperty("color", mitk::ColorProperty::New(b2, b1, b0)); GetDefaultDataStorage()->Add(node, m_SelectedDiffusionNodes.front()); shellCount++; } // add origin to datastorage mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(originSet); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); name += "_Origin"; node->SetName(name.toStdString().c_str()); node->SetProperty("pointsize", mitk::FloatProperty::New((float)maxSize/50)); node->SetProperty("color", mitk::ColorProperty::New(1,1,1)); GetDefaultDataStorage()->Add(node, m_SelectedDiffusionNodes.front()); } void QmitkPreprocessingView::DoReduceGradientDirections() { if (m_DiffusionImage.IsNull()) return; typedef mitk::DiffusionImage DiffusionImageType; typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter FilterType; typedef DiffusionImageType::BValueMap BValueMap; // GetShellSelection from GUI BValueMap shellSlectionMap; BValueMap originalShellMap = m_DiffusionImage->GetB_ValueMap(); std::vector newNumGradientDirections; int shellCounter = 0; foreach(QCheckBox * box , m_ReduceGradientCheckboxes) { if(box->isChecked()) { double BValue = (box->text().split(' ')).at(0).toDouble(); shellSlectionMap[BValue] = originalShellMap[BValue]; newNumGradientDirections.push_back(m_ReduceGradientSpinboxes.at(shellCounter)->value()); } shellCounter++; } if (newNumGradientDirections.empty()) return; GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetDirections(); FilterType::Pointer filter = FilterType::New(); filter->SetInput(m_DiffusionImage->GetVectorImage()); filter->SetOriginalGradientDirections(gradientContainer); filter->SetNumGradientDirections(newNumGradientDirections); filter->SetOriginalBValueMap(originalShellMap); filter->SetShellSelectionBValueMap(shellSlectionMap); filter->Update(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(m_DiffusionImage->GetB_Value()); image->SetDirections(filter->GetGradientDirections()); image->SetMeasurementFrame(m_DiffusionImage->GetMeasurementFrame()); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); foreach(QSpinBox* box, m_ReduceGradientSpinboxes) { name += "_"; name += QString::number(box->value()); } imageNode->SetName(name.toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); } void QmitkPreprocessingView::MergeDwis() { typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; if (m_SelectedDiffusionNodes.size()<2) return; typedef itk::VectorImage DwiImageType; typedef DwiImageType::PixelType DwiPixelType; typedef DwiImageType::RegionType DwiRegionType; typedef std::vector< DwiImageType::Pointer > DwiImageContainerType; typedef std::vector< GradientContainerType::Pointer > GradientListContainerType; DwiImageContainerType imageContainer; GradientListContainerType gradientListContainer; std::vector< double > bValueContainer; QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); for (int i=0; i* >( m_SelectedDiffusionNodes.at(i)->GetData() ); if ( dwi.IsNotNull() ) { imageContainer.push_back(dwi->GetVectorImage()); gradientListContainer.push_back(dwi->GetDirections()); bValueContainer.push_back(dwi->GetB_Value()); if (i>0) { name += "+"; name += m_SelectedDiffusionNodes.at(i)->GetName().c_str(); } } } typedef itk::MergeDiffusionImagesFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetImageVolumes(imageContainer); filter->SetGradientLists(gradientListContainer); filter->SetBValues(bValueContainer); filter->Update(); vnl_matrix_fixed< double, 3, 3 > mf; mf.set_identity(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(filter->GetB_Value()); image->SetDirections(filter->GetOutputGradients()); image->SetMeasurementFrame(mf); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); imageNode->SetName(name.toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); } void QmitkPreprocessingView::ExtractB0() { typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; int nrFiles = m_SelectedDiffusionNodes.size(); if (!nrFiles) return; // call the extraction withou averaging if the check-box is checked if( this->m_Controls->m_CheckExtractAll->isChecked() ) { DoExtractBOWithoutAveraging(); return; } mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes.begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes.end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { DiffusionImageType* vols = static_cast( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); filter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( filter->GetOutput() ); mitkImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_B0")); GetDefaultDataStorage()->Add(node); ++itemiter; } } void QmitkPreprocessingView::DoExtractBOWithoutAveraging() { // typedefs typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; typedef itk::B0ImageExtractionToSeparateImageFilter< short, short> FilterType; // check number of selected objects, return if empty int nrFiles = m_SelectedDiffusionNodes.size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes.begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes.end() ); while ( itemiter != itemiterend ) // for all items { DiffusionImageType* vols = static_cast( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); filter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( filter->GetOutput() ); mitkImage->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_B0_ALL")); GetDefaultDataStorage()->Add(node); ++itemiter; } } void QmitkPreprocessingView::AverageGradients() { int nrFiles = m_SelectedDiffusionNodes.size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes.begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes.end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); vols->AverageRedundantGradients(m_Controls->m_Blur->value()); ++itemiter; } }