diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp index f9d2b4f574..dead704265 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.cpp @@ -1,396 +1,336 @@ /*=================================================================== 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_MultiShellDirectionalKurtosisFitImageFilter_cpp_ #define _itk_MultiShellDirectionalKurtosisFitImageFilter_cpp_ #endif #define _USE_MATH_DEFINES #include "itkMultiShellRadialAdcKurtosisImageFilter.h" #include #include #include "mitkDiffusionFunctionCollection.h" -#include "vnl/vnl_least_squares_function.h" -#include "vnl/algo/vnl_levenberg_marquardt.h" - -namespace itk -{ -/** baseclass for IVIM fitting algorithms */ -struct lestSquaresFunction: public vnl_least_squares_function +namespace itk { - void set_measurements(const vnl_vector& x) - { - measurements.set_size(x.size()); - measurements.copy_in(x.data_block()); - } - - void set_bvalues(const vnl_vector& x) - { - b.set_size(x.size()); - b.copy_in(x.data_block()); - } - - void set_reference_measuremnt(const double & x) - { - reference_measurement = x; - } - - double reference_measurement; - vnl_vector measurements; - vnl_vector b; - - int N; - - lestSquaresFunction(unsigned int number_of_measurements) : - vnl_least_squares_function(2, number_of_measurements, no_gradient) - { - N = get_number_of_residuals(); - } - - void f(const vnl_vector& x, vnl_vector& fx) { - - const double & D = x[0]; - const double & K = x[1]; - - for(int s=0; s MultiShellRadialAdcKurtosisImageFilter ::MultiShellRadialAdcKurtosisImageFilter() { this->SetNumberOfRequiredInputs( 1 ); - this->SetNumberOfThreads(1); } template void MultiShellRadialAdcKurtosisImageFilter ::BeforeThreadedGenerateData() { - // test whether BvalueMap contains all necessary information if(m_B_ValueMap.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_B_Value * gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; m_B_ValueMap[bValueKey].push_back(gdcit.Index()); } } //# BValueMap contains no bZero --> itkException if(m_B_ValueMap.find(0.0) == m_B_ValueMap.end()) { MITK_INFO << "No ReferenceSignal (BZeroImages) found!"; itkExceptionMacro(<< "No ReferenceSignal (BZeroImages) found!"); } // [allDirectionsContainer] Gradient DirectionContainer containing all unique directions m_allDirectionsIndicies = mitk::gradients::GetAllUniqueDirections(m_B_ValueMap, m_OriginalGradientDirections); // [sizeAllDirections] size of GradientContainer cointaining all unique directions m_allDirectionsSize = m_allDirectionsIndicies.size(); m_TargetGradientDirections = mitk::gradients::CreateNormalizedUniqueGradientDirectionContainer(m_B_ValueMap,m_OriginalGradientDirections); m_ShellInterpolationMatrixVector.reserve(m_B_ValueMap.size()-1); - m_bValueVector = vnl_vector(2); + m_bValueVector = vnl_vector(m_B_ValueMap.size()-1); // for each shell BValueMap::const_iterator it = m_B_ValueMap.begin(); it++; //skip bZeroIndices unsigned int shellIndex = 0; - - for(;it != m_B_ValueMap.end();it++) + for(;it != m_B_ValueMap.end();++it) { //- calculate maxShOrder const IndicesVector currentShell = it->second; unsigned int SHMaxOrder = 12; while( ((SHMaxOrder+1)*(SHMaxOrder+2)/2) > currentShell.size() && ((SHMaxOrder+1)*(SHMaxOrder+2)/2) >= 4 ) - { - SHMaxOrder -= 2 ; - } + SHMaxOrder -= 2 ; + //- get TragetSHBasis using allDirectionsContainer vnl_matrix sphericalCoordinates; sphericalCoordinates = mitk::gradients::ComputeSphericalFromCartesian(m_allDirectionsIndicies, m_OriginalGradientDirections); vnl_matrix TargetSHBasis = mitk::gradients::ComputeSphericalHarmonicsBasis(sphericalCoordinates, SHMaxOrder); - //MITK_INFO << "TargetSHBasis " << TargetSHBasis.rows() << " x " << TargetSHBasis.cols(); //- get ShellSHBasis using currentShellDirections sphericalCoordinates = mitk::gradients::ComputeSphericalFromCartesian(currentShell, m_OriginalGradientDirections); vnl_matrix ShellSHBasis = mitk::gradients::ComputeSphericalHarmonicsBasis(sphericalCoordinates, SHMaxOrder); - //MITK_INFO << "ShellSHBasis " << ShellSHBasis.rows() << " x " << ShellSHBasis.cols(); //- calculate interpolationSHBasis [TargetSHBasis * ShellSHBasis^-1] vnl_matrix_inverse invShellSHBasis(ShellSHBasis); vnl_matrix shellInterpolationMatrix = TargetSHBasis * invShellSHBasis.pinverse(); - //MITK_INFO << "shellInterpolationMatrix " << shellInterpolationMatrix.rows() << " x " << shellInterpolationMatrix.cols(); //- save interpolationSHBasis m_ShellInterpolationMatrixVector.push_back(shellInterpolationMatrix); m_bValueVector.put(shellIndex,it->first); ++shellIndex; } - m_WeightsVector.reserve(m_B_ValueMap.size()-1); BValueMap::const_iterator itt = m_B_ValueMap.begin(); itt++; // skip ReferenceImages //- calculate Weights [Weigthing = shell_size / max_shell_size] unsigned int maxShellSize = 0; for(;itt != m_B_ValueMap.end(); itt++){ if(itt->second.size() > maxShellSize) maxShellSize = itt->second.size(); } itt = m_B_ValueMap.begin(); itt++; // skip ReferenceImages for(;itt != m_B_ValueMap.end(); itt++){ m_WeightsVector.push_back(itt->second.size() / (double)maxShellSize); MITK_INFO << m_WeightsVector.back(); } - // initialize output image typename OutputImageType::Pointer outImage = static_cast(ProcessObject::GetOutput(0)); 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( 1+m_allDirectionsSize ); // size of 1(bzeroValue) + AllDirectionsContainer outImage->Allocate(); BValueMap::iterator ittt = m_B_ValueMap.begin(); ittt++; // skip bZeroImages corresponding to 0-bValue m_TargetB_Value = 0; while(ittt!=m_B_ValueMap.end()) { m_TargetB_Value += ittt->first; ittt++; } m_TargetB_Value /= (double)(m_B_ValueMap.size()-1); MITK_INFO << "Input:" << std::endl << std::endl << " GradientDirections: " << m_OriginalGradientDirections->Size() << std::endl << " Shells: " << (m_B_ValueMap.size() - 1) << std::endl << " ReferenceImages: " << m_B_ValueMap[0.0].size() << std::endl; MITK_INFO << "Output:" << std::endl << std::endl << " OutImageVectorLength: " << outImage->GetVectorLength() << std::endl << " TargetDirections: " << m_allDirectionsSize << std::endl << " TargetBValue: " << m_TargetB_Value << std::endl << std::endl; } template void MultiShellRadialAdcKurtosisImageFilter ::ThreadedGenerateData(const OutputImageRegionType &outputRegionForThread, ThreadIdType /*threadId*/) { - // 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(); // calculate target bZero-Value [b0_t] const IndicesVector BZeroIndices = m_B_ValueMap[0.0]; double BZeroAverage = 0.0; unsigned int numberOfShells = m_B_ValueMap.size()-1; // create empty nxm SignalMatrix containing n->signals/directions (in case of interpolation ~ sizeAllDirections otherwise the size of any shell) for m->shells vnl_matrix SignalMatrix(m_allDirectionsSize, numberOfShells); // create nx1 targetSignalVector vnl_vector SignalVector(m_allDirectionsSize); OutputPixelType out; InputPixelType b; BValueMap::const_iterator shellIterator; vnl_matrix lsfCoeffs; vnl_vector InterpVector; // ** walking over each Voxel while(!iit.IsAtEnd()) { b = iit.Get(); BZeroAverage=0.0; for(unsigned int i = 0 ; i < BZeroIndices.size(); i++){ BZeroAverage += b[BZeroIndices[i]]; } BZeroAverage /= (double)BZeroIndices.size(); out = oit.Get(); out.Fill(0.0); out.SetElement(0,BZeroAverage); shellIterator = m_B_ValueMap.begin(); shellIterator++; //skip bZeroImages unsigned int shellIndex = 0; while(shellIterator != m_B_ValueMap.end()) { // reset Data SignalVector.fill(0.0); // - get the RawSignal const IndicesVector & currentShell = shellIterator->second; InterpVector.set_size(currentShell.size()); // - get raw Signal for currente shell for(unsigned int i = 0 ; i < currentShell.size(); i++) InterpVector.put(i,b[currentShell[i]]); //- interpolate the Signal using corresponding interpolationSHBasis SignalVector = m_ShellInterpolationMatrixVector.at(shellIndex) * InterpVector; SignalMatrix.set_column(shellIndex, SignalVector); shellIterator++; shellIndex++; } - lsfCoeffs.set_size(m_allDirectionsSize , 2); - calculateCoeffs(lsfCoeffs,SignalMatrix, m_bValueVector, BZeroAverage); - calculateSignalFromLsfCoeffs(SignalVector,lsfCoeffs,m_TargetB_Value,BZeroAverage); + lsfCoeffs.set_size(m_allDirectionsSize , 2 /*Number of unknowns [ADC AKC]*/); + calculateCoeffs(lsfCoeffs, SignalMatrix, BZeroAverage); + calculateSignalFromLsfCoeffs(SignalVector,lsfCoeffs, BZeroAverage, m_TargetB_Value); for(unsigned int i = 1 ; i < out.Size(); i ++) out.SetElement(i,SignalVector.get(i-1)); oit.Set(out); ++oit; ++iit; } } template void MultiShellRadialAdcKurtosisImageFilter ::logarithm( vnl_vector & vec) { for(unsigned int i = 0; i < vec.size(); i++){ if(vec[i] <= 0.0) vec[i] = 0.00001; vec[i] = std::log( vec[i] ); } } - template void MultiShellRadialAdcKurtosisImageFilter -::calculateCoeffs(vnl_matrix &lsfCoeffs, const vnl_matrix & SignalMatrix, const vnl_vector & bValueVector, const double & reference_b_value) +::calculateCoeffs(vnl_matrix &lsfCoeffs, const vnl_matrix & SignalMatrix, const double & S0) { vnl_vector initalGuess(2); // initialize Least Squres Function - lestSquaresFunction model(bValueVector.size()); + // SignalMatrix.cols() defines the number of shells/measurement points + lestSquaresFunction model(SignalMatrix.cols()); + // set BValue Vector e.g.: [1000, 2000, 3000] <- shell b Values + model.set_bvalues(m_bValueVector); + // initialize Levenberg Marquardt vnl_levenberg_marquardt minimizer(model); + minimizer.set_max_function_evals(10000); // Iterations minimizer.set_f_tolerance(1e-10); // Function tolerance - minimizer.set_epsilon_function(1e-10); // Epsilon // for each Direction calculate LSF Coeffs ADC & AKC for(unsigned int i = 0 ; i < SignalMatrix.rows(); i++) { - // set Signal Vector + // set voxel signal vector model.set_measurements(SignalMatrix.get_row(i)); - // set BValue Vector e.g.: [1000, 2000, 3000] - model.set_bvalues(bValueVector); - // set BZero Value - model.set_reference_measuremnt(reference_b_value); - + // set voxel reference signal + model.set_reference_measurement(S0); // Start Vector [ADC, AKC] - initalGuess.put(0, 0.f); - initalGuess.put(1, 0.f); + initalGuess.put(0, 0.f); // ADC + initalGuess.put(1, 0.8f); // AKC // start Levenberg-Marquardt bool status = minimizer.minimize_without_gradient(initalGuess); if(!status) { MITK_INFO<< "______________________________"; MITK_INFO<< "Minimizer f Error: " << minimizer.get_f_tolerance(); MITK_INFO<< "Minimizer end Error: " << minimizer.get_end_error(); MITK_INFO<< "Minimizer Evaluations: " << minimizer.get_num_evaluations(); MITK_INFO<< "Minimizer Iterations: " << minimizer.get_num_iterations(); MITK_INFO<< "______________________________"; - }else{ - std::cout << SignalMatrix.get_row(i)[0] << ";" // VALVEC - << SignalMatrix.get_row(i)[1] << ";" - << SignalMatrix.get_row(i)[2] << ";" - << reference_b_value << ";"//REFVAL - << initalGuess[0] << ";" // ADC - << initalGuess[1] << std::endl; // AKC - } + }/*else{ //OUTPUT FOR EVALUATION + + std::cout << std::scientific << std::setprecision(5) + << initalGuess[0] << ";" // fitted ADC + << initalGuess[1] << ";" // fitted AKC + << S0 << ";" // S0 value + << minimizer.get_end_error() << ";"; // End error + for(unsigned int j = 0; j < SignalMatrix.get_row(i).size(); j++ ) + std::cout << std::scientific << std::setprecision(5) + << SignalMatrix.get_row(i)[j] << ";"; // S_n Values corresponding to shell 1 to shell n + std::cout << std::endl; + }*/ // Set Coeffs for each gradient direction lsfCoeffs.set_row(i, initalGuess); } } - template void MultiShellRadialAdcKurtosisImageFilter -::calculateSignalFromLsfCoeffs(vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double &bValue, const double & referenceSignal) +::calculateSignalFromLsfCoeffs(vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & S0, const double &b /*target bValue*/) { // For each Direction for(unsigned int i = 0 ; i < lsfCoeffs.rows();i++){ const double & D = lsfCoeffs(i,0); const double & K = lsfCoeffs(i,1); - - vec[i] = referenceSignal * exp( -bValue * D + 1./6.* bValue * bValue * D * D * K); - - std::cout << "OUTVAL = " << i << std::endl; + vec[i] = S0 * exp( -b * D + 1./6.* b * b * D * D * K); } } - } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h index 50192e1e19..3326cb1cd2 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellRadialAdcKurtosisImageFilter.h @@ -1,120 +1,167 @@ /*=================================================================== 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_MultiShellRadialAdcKurtosisImageFilter_h_ #define _itk_MultiShellRadialAdcKurtosisImageFilter_h_ #include #include #include +#include "vnl/vnl_least_squares_function.h" +#include "vnl/algo/vnl_levenberg_marquardt.h" namespace itk { /** * \brief Select subset of the input vectors equally distributed over the sphere using an iterative electrostatic repulsion strategy. */ template class MultiShellRadialAdcKurtosisImageFilter : public ImageToImageFilter, itk::VectorImage > { + private: + struct lestSquaresFunction: public vnl_least_squares_function + { + + void set_measurements(const vnl_vector& x) + { + measurements.set_size(x.size()); + measurements.copy_in(x.data_block()); + } + + void set_bvalues(const vnl_vector& x) + { + bValueVector.set_size(x.size()); + bValueVector.copy_in(x.data_block()); + } + + void set_reference_measurement(const double & x) + { + S0 = x; + } + + vnl_vector measurements; + vnl_vector bValueVector; + double S0; + int N; + + lestSquaresFunction(unsigned int number_of_measurements) : + vnl_least_squares_function(2 /*number of unknowns [ADC AKC]*/, number_of_measurements, no_gradient) + { + N = get_number_of_residuals(); + } + + void f(const vnl_vector& x, vnl_vector& fx) { + + const double & D = x[0]; + const double & K = x[1]; + const vnl_vector & b = bValueVector; + + for(int s=0; s 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: MultiShellRadialAdcKurtosisImageFilter(); ~MultiShellRadialAdcKurtosisImageFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads ); void logarithm( vnl_vector & vec); - void calculateCoeffs( vnl_matrix & lsfCoeffs, const vnl_matrix & SignalMatrix, const vnl_vector & bValueVector, const double & reference_b_value); - void calculateSignalFromLsfCoeffs( vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & bValue, const double & referenceSignal); - - + void calculateCoeffs( vnl_matrix & lsfCoeffs, const vnl_matrix & SignalMatrix, const double & S0); + void calculateSignalFromLsfCoeffs( vnl_vector & vec, const vnl_matrix & lsfCoeffs, const double & S0, const double & b); 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; vnl_vector m_bValueVector; std::vector m_bZeroIndicesSplitVectors; IndicesVector m_allDirectionsIndicies; unsigned int m_allDirectionsSize; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkMultiShellRadialAdcKurtosisImageFilter.cpp" #endif #endif 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 86a82a93ad..d1483b9f67 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp @@ -1,730 +1,771 @@ /*=================================================================== 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 +#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 #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()) ); + connect( (QObject*)(m_Controls->m_AdcSignalFit), SIGNAL(clicked()), this, SLOT(DoADCFit()) ); + connect( (QObject*)(m_Controls->m_AkcSignalFit), SIGNAL(clicked()), this, SLOT(DoAKCFit()) ); connect( (QObject*)(m_Controls->m_B_ValueMap_Rounder_SpinBox), SIGNAL(valueChanged(int)), this, SLOT(UpdateDwiBValueMapRounder(int))); connect( (QObject*)(m_Controls->m_CreateLengthCorrectedDwi), SIGNAL(clicked()), this, SLOT(DoLengthCorrection()) ); connect( (QObject*)(m_Controls->m_CalcAdcButton), SIGNAL(clicked()), this, SLOT(DoAdcCalculation()) ); } } void QmitkPreprocessingView::DoLengthCorrection() { if (m_DiffusionImage.IsNull()) return; typedef mitk::DiffusionImage DiffusionImageType; typedef itk::DwiGradientLengthCorrectionFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetRoundingValue( m_Controls->m_B_ValueMap_Rounder_SpinBox->value()); filter->SetReferenceBValue(m_DiffusionImage->GetB_Value()); filter->SetReferenceGradientDirectionContainer(m_DiffusionImage->GetDirections()); filter->Update(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( m_DiffusionImage->GetVectorImage()); image->SetB_Value( filter->GetNewBValue() ); image->SetDirections( filter->GetOutputGradientDirectionContainer()); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); imageNode->SetName((name+"_rounded").toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); } void QmitkPreprocessingView::UpdateDwiBValueMapRounder(int i) { if (m_DiffusionImage.IsNull()) return; m_DiffusionImage->UpdateBValueMap(); UpdateBValueTableWidget(i); } -void QmitkPreprocessingView::DoAdcAverage() +void QmitkPreprocessingView::DoBiExpFit(){} + +void QmitkPreprocessingView::DoAKCFit() +{ + typedef mitk::DiffusionImage DiffusionImageType; + typedef itk::MultiShellRadialAdcKurtosisImageFilter FilterType; + typedef DiffusionImageType::BValueMap BValueMap; + + for (int i=0; i(m_SelectedDiffusionNodes.at(i)->GetData()); + BValueMap originalShellMap = inImage->GetB_ValueMap(); + + GradientDirectionContainerType::Pointer gradientContainer = inImage->GetDirections(); + + FilterType::Pointer filter = FilterType::New(); + filter->SetInput(inImage->GetVectorImage()); + filter->SetOriginalGradientDirections(gradientContainer); + filter->SetOriginalBValueMap(originalShellMap); + filter->SetB_Value(inImage->GetB_Value()); + filter->Update(); + + DiffusionImageType::Pointer outImage = DiffusionImageType::New(); + outImage->SetVectorImage( filter->GetOutput() ); + outImage->SetB_Value( filter->GetTargetB_Value() ); + outImage->SetDirections( filter->GetTargetGradientDirections() ); + outImage->InitializeFromVectorImage(); + + mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); + imageNode->SetData( outImage ); + QString name = m_SelectedDiffusionNodes.at(i)->GetName().c_str(); + + imageNode->SetName((name+"_averaged").toStdString().c_str()); + GetDefaultDataStorage()->Add(imageNode); + } +} + + +void QmitkPreprocessingView::DoADCFit() { typedef mitk::DiffusionImage DiffusionImageType; typedef itk::MultiShellAdcAverageReconstructionImageFilter FilterType; typedef DiffusionImageType::BValueMap BValueMap; for (int i=0; i(m_SelectedDiffusionNodes.at(i)->GetData()); BValueMap originalShellMap = inImage->GetB_ValueMap(); GradientDirectionContainerType::Pointer gradientContainer = inImage->GetDirections(); FilterType::Pointer filter = FilterType::New(); filter->SetInput(inImage->GetVectorImage()); filter->SetOriginalGradientDirections(gradientContainer); filter->SetOriginalBValueMap(originalShellMap); filter->SetB_Value(inImage->GetB_Value()); filter->Update(); DiffusionImageType::Pointer outImage = DiffusionImageType::New(); outImage->SetVectorImage( filter->GetOutput() ); outImage->SetB_Value( filter->GetTargetB_Value() ); //image->SetOriginalDirections( filter->GetTargetGradientDirections() ); outImage->SetDirections( filter->GetTargetGradientDirections() ); //image->SetMeasurementFrame( m_DiffusionImage->GetMeasurementFrame() ); outImage->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( outImage ); QString name = m_SelectedDiffusionNodes.at(i)->GetName().c_str(); imageNode->SetName((name+"_averaged").toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); } } void QmitkPreprocessingView::DoAdcCalculation() { if (m_DiffusionImage.IsNull()) return; typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; typedef itk::AdcImageFilter< DiffusionPixelType, double > FilterType; for (int i=0; i(m_SelectedDiffusionNodes.at(i)->GetData()); FilterType::Pointer filter = FilterType::New(); filter->SetInput(inImage->GetVectorImage()); filter->SetGradientDirections(inImage->GetDirections()); filter->SetB_value(inImage->GetB_Value()); filter->Update(); mitk::Image::Pointer image = mitk::Image::New(); image->InitializeByItk( filter->GetOutput() ); image->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); QString name = m_SelectedDiffusionNodes.at(i)->GetName().c_str(); imageNode->SetName((name+"_ADC").toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); } } void QmitkPreprocessingView::UpdateBValueTableWidget(int i) { 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 (m_DiffusionImage.IsNull()) { m_Controls->m_B_ValueMap_TableWidget->clear(); m_Controls->m_B_ValueMap_TableWidget->setRowCount(1); QStringList headerList; headerList << "b-Value" << "Number of gradients"; m_Controls->m_B_ValueMap_TableWidget->setHorizontalHeaderLabels(headerList); m_Controls->m_B_ValueMap_TableWidget->setItem(0,0,new QTableWidgetItem("-")); m_Controls->m_B_ValueMap_TableWidget->setItem(0,1,new QTableWidgetItem("-")); }else{ typedef mitk::DiffusionImage::BValueMap BValueMap; typedef mitk::DiffusionImage::BValueMap::iterator BValueMapIterator; BValueMapIterator it; BValueMap roundedBValueMap; GradientDirectionContainerType::ConstIterator gdcit; for( gdcit = m_DiffusionImage->GetDirections()->Begin(); gdcit != m_DiffusionImage->GetDirections()->End(); ++gdcit) { float currentBvalue = std::floor(m_DiffusionImage->GetB_Value(gdcit.Index())); unsigned int rounded = int((currentBvalue + 0.5 * i)/i)*i; roundedBValueMap[rounded].push_back(gdcit.Index()); } m_Controls->m_B_ValueMap_TableWidget->clear(); m_Controls->m_B_ValueMap_TableWidget->setRowCount(roundedBValueMap.size() ); QStringList headerList; headerList << "b-Value" << "Number of gradients"; m_Controls->m_B_ValueMap_TableWidget->setHorizontalHeaderLabels(headerList); QCheckBox* checkBox; QSpinBox* spinBox; int i = 0 ; for(it = roundedBValueMap.begin() ;it != roundedBValueMap.end(); it++) { m_Controls->m_B_ValueMap_TableWidget->setItem(i,0,new QTableWidgetItem(QString::number(it->first))); m_Controls->m_B_ValueMap_TableWidget->setItem(i,1,new QTableWidgetItem(QString::number(it->second.size()))); // Reduce Gradients GUI adaption if(it->first != 0 && roundedBValueMap.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++; } } } 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); m_Controls->m_B_ValueMap_Rounder_SpinBox->setEnabled(foundDwiVolume); - m_Controls->m_AdcAverage->setEnabled(foundDwiVolume); + m_Controls->m_AdcSignalFit->setEnabled(foundDwiVolume); + m_Controls->m_AkcSignalFit->setEnabled(foundDwiVolume); m_Controls->m_CreateLengthCorrectedDwi->setEnabled(foundDwiVolume); m_Controls->m_CalcAdcButton->setEnabled(foundDwiVolume); UpdateBValueTableWidget(m_Controls->m_B_ValueMap_Rounder_SpinBox->value()); 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); } } 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_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; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h index 5bfeb97b89..a44e6018a7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h @@ -1,115 +1,117 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _QMITKPREPROCESSINGVIEW_H_INCLUDED #define _QMITKPREPROCESSINGVIEW_H_INCLUDED #include #include #include "ui_QmitkPreprocessingViewControls.h" #include "mitkDiffusionImage.h" typedef short DiffusionPixelType; struct PrpSelListener; /*! * \ingroup org_mitk_gui_qt_preprocessing_internal * * \brief QmitkPreprocessingView * * Document your class here. * * \sa QmitkFunctionality */ class QmitkPreprocessingView : public QmitkFunctionality { friend struct PrpSelListener; // this is needed for all Qt objects that should have a MOC object (everything that derives from QObject) Q_OBJECT public: static const std::string VIEW_ID; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; QmitkPreprocessingView(); QmitkPreprocessingView(const QmitkPreprocessingView& other); virtual ~QmitkPreprocessingView(); virtual void CreateQtPartControl(QWidget *parent); /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); /// \brief Called when the functionality is activated virtual void Activated(); virtual void Deactivated(); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); static const int nrconvkernels; protected slots: void AverageGradients(); void ExtractB0(); void MergeDwis(); void DoApplyMesurementFrame(); void DoReduceGradientDirections(); void DoShowGradientDirections(); void DoHalfSphereGradientDirections(); - void DoAdcAverage(); + void DoADCFit(); + void DoAKCFit(); + void DoBiExpFit(); void UpdateDwiBValueMapRounder(int i); void DoLengthCorrection(); void DoAdcCalculation(); protected: /** Called by ExtractB0 if check-box activated, extracts all b0 images without averaging */ void DoExtractBOWithoutAveraging(); void UpdateBValueTableWidget(int i); /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); Ui::QmitkPreprocessingViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; void SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name); mitk::DiffusionImage::Pointer m_DiffusionImage; std::vector< mitk::DataNode::Pointer > m_SelectedDiffusionNodes; QList m_ReduceGradientCheckboxes; QList m_ReduceGradientSpinboxes; }; #endif // _QMITKPREPROCESSINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui index af173bb62d..a3c921ef5e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui @@ -1,648 +1,662 @@ QmitkPreprocessingViewControls 0 0 892 1306 0 0 false QmitkPreprocessingViewControls true Please Select Input Data Raw DWI: <html><head/><body><p><span style=" color:#ff0000;">mandatory</span></p></body></html> true 0 0 Info - - + + false + + Generate pointset displaying the gradient vectors (applied measurement frame). + + + + + + + - Round b-value + Show gradients 0 0 Qt::ScrollBarAsNeeded Qt::ScrollBarAlwaysOff true 100 true false true b-Value Number of gradients - - + + false - - Generate pointset displaying the gradient vectors (applied measurement frame). - - - - - - - - Show gradients + Round b-value Qt::Horizontal Qt::Horizontal Qt::Horizontal - - - - false - - - <html><head/><body><p>Create a new Diffusion Weighted Image (DWI) using an ADC average over all q-shells.</p></body></html> - - - ADC average (multi-shell) - - - <html><head/><body><p>Define the sampling frame the b-Values are rounded with.</p></body></html> Sampling frame: false <html><head/><body><p>Round b-values to nearest multiple of this value (click &quot;Round b-value&quot; to create new image with these values).</p></body></html> QAbstractSpinBox::CorrectToNearestValue 1 10000 10 + + + + + + false + + + <html><head/><body><p>Create a new Diffusion Weighted Image (DWI) using an ADC average over all q-shells.</p></body></html> + + + ADC Signal Fit + + + + + + + false + + + AKC Signal Fit + + + + + Non diffusion weighted image 0 30 Average and extract all images that were acquired without diffusion weighting. true false If multiple baseline acquisitions are present, the default behaviour is to output an averaged image. Extract B0 Create a 3D+t data set containing all b0 images as timesteps Extract all B0 images without averaging false If multiple baseline acquisitions are present, the default behaviour is to output an averaged image. Calculate ADC map Modify DWI QFrame::NoFrame QFrame::Raised 0 Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. 6 2.000000000000000 0.000100000000000 0.001000000000000 Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Merge radius false Multiple acquistions of one gradient direction can be averaged. Due to rounding errors, similar gradients often differ in the last decimal positions. The Merge radius allows to average them by taking all directions within a certain radius into account. Average redundant gradients Qt::Horizontal QFrame::NoFrame QFrame::Raised 0 9 Specify desired number of gradients per shell: Qt::Horizontal 40 20 false Retain only the specified number of gradient directions and according image volumes. The retained directions are spread equally over the half sphere using an iterative energy repulsion strategy. Reduce number of gradients Qt::Horizontal false Sometimes the gradient directions are not located on one half sphere. Mirror gradients to half sphere Merge selected images false Merges selected DWIs of same dimension. If several b-values are present, the resulting image will contain multiple b-shells. Merge selected DWIs 0 0 Measurment frame Qt::Horizontal 40 20 false 0 0 0 0 IBeamCursor true Qt::ScrollBarAlwaysOff Qt::ScrollBarAlwaysOff true false false true true 0 false true true New Row New Row New Row New Column New Column New Column false Apply new measurement frame Qt::Vertical 20 40