diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp index b166ffabd5..21ee40e231 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.cpp @@ -1,324 +1,324 @@ /*=================================================================== 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 "mitkDiffusionFunctionCollection.h" namespace itk { template MultiShellAdcAverageReconstructionImageFilter ::MultiShellAdcAverageReconstructionImageFilter() { this->SetNumberOfRequiredInputs( 1 ); // this->SetNumberOfThreads(1); } template void MultiShellAdcAverageReconstructionImageFilter ::BeforeThreadedGenerateData() { // test whether BvalueMap contains all necessary information - if(m_BValueMap.size() == 0) + 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_BValue * gdcit.Value().two_norm() * gdcit.Value().two_norm())+7.5)/10)*10; - m_BValueMap[bValueKey].push_back(gdcit.Index()); + 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_BValueMap.find(0.0) == m_BValueMap.end()) + 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_BValueMap, m_OriginalGradientDirections); + 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_BValueMap,m_OriginalGradientDirections); + m_TargetGradientDirections = mitk::gradients::CreateNormalizedUniqueGradientDirectionContainer(m_B_ValueMap,m_OriginalGradientDirections); - for(BValueMap::const_iterator it = m_BValueMap.begin();it != m_BValueMap.end(); it++) + for(BValueMap::const_iterator it = m_B_ValueMap.begin();it != m_B_ValueMap.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"); } } - m_ShellInterpolationMatrixVector.reserve(m_BValueMap.size()-1); + m_ShellInterpolationMatrixVector.reserve(m_B_ValueMap.size()-1); // for each shell - BValueMap::const_iterator it = m_BValueMap.begin(); + BValueMap::const_iterator it = m_B_ValueMap.begin(); it++; //skip bZeroIndices - for(;it != m_BValueMap.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 -= 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_WeightsVector.reserve(m_BValueMap.size()-1); - BValueMap::const_iterator itt = m_BValueMap.begin(); + 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_BValueMap.end(); itt++){ + for(;itt != m_B_ValueMap.end(); itt++){ if(itt->second.size() > maxShellSize) maxShellSize = itt->second.size(); } - itt = m_BValueMap.begin(); + itt = m_B_ValueMap.begin(); itt++; // skip ReferenceImages - for(;itt != m_BValueMap.end(); itt++){ + 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 = 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( 1+m_allDirectionsSize ); // size of 1(bzeroValue) + AllDirectionsContainer outImage->Allocate(); - BValueMap::iterator ittt = m_BValueMap.begin(); + BValueMap::iterator ittt = m_B_ValueMap.begin(); ittt++; // skip bZeroImages corresponding to 0-bValue - m_TargetBValue = 0; - while(ittt!=m_BValueMap.end()) + m_TargetB_Value = 0; + while(ittt!=m_B_ValueMap.end()) { - m_TargetBValue += ittt->first; + m_TargetB_Value += ittt->first; ittt++; } - m_TargetBValue /= (double)(m_BValueMap.size()-1); + m_TargetB_Value /= (double)(m_B_ValueMap.size()-1); MITK_INFO << "Input:" << std::endl << std::endl << " GradientDirections: " << m_OriginalGradientDirections->Size() << std::endl - << " Shells: " << (m_BValueMap.size() - 1) << std::endl - << " ReferenceImages: " << m_BValueMap.at(0.0).size() << std::endl; + << " Shells: " << (m_B_ValueMap.size() - 1) << std::endl + << " ReferenceImages: " << m_B_ValueMap.at(0.0).size() << std::endl; MITK_INFO << "Output:" << std::endl << std::endl << " OutImageVectorLength: " << outImage->GetVectorLength() << std::endl << " TargetDirections: " << m_allDirectionsSize << std::endl - << " TargetBValue: " << m_TargetBValue << std::endl + << " TargetBValue: " << m_TargetB_Value << std::endl << std::endl; } template void MultiShellAdcAverageReconstructionImageFilter ::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_BValueMap[0.0]; + const IndicesVector BZeroIndices = m_B_ValueMap[0.0]; double BZeroAverage = 0.0; // 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, m_BValueMap.size()-1); + vnl_matrix SignalMatrix(m_allDirectionsSize, m_B_ValueMap.size()-1); // create nx1 targetSignalVector vnl_vector SignalVector(m_allDirectionsSize); // ** walking over each Voxel while(!iit.IsAtEnd()) { InputPixelType b = iit.Get(); BZeroAverage=0.0; for(unsigned int i = 0 ; i < BZeroIndices.size(); i++){ //MITK_INFO << "BValues("<second; vnl_vector InterpVector(currentShell.size()); // - get raw Signal for currente shell for(unsigned int i = 0 ; i < currentShell.size(); i++) InterpVector.put(i,b[currentShell[i]]); //MITK_INFO << "InputVector(RAWSIGNAL:"<< shellIndex <<")"; //MITK_INFO << InterpVector ; //- normalization of the raw Signal S_S0Normalization(InterpVector, BZeroAverage); //MITK_INFO << "InputVector(NORMALIZED:"<< shellIndex <<")"; //MITK_INFO << InterpVector ; //- interpolate the Signal if necessary using corresponding interpolationSHBasis SignalVector = m_ShellInterpolationMatrixVector.at(shellIndex) * InterpVector; //MITK_INFO << "SignalVector(NORMALIZED:"<< shellIndex <<")"; //MITK_INFO << SignalVector ; // - ADC calculation for the signalVector calculateAdcFromSignal(SignalVector, shellIterator->first); //MITK_INFO << "SignalVector(SIGNaL->ADC:"<< shellIndex <<")"; //MITK_INFO << SignalVector ; //- weight the signal SignalVector *= m_WeightsVector.at(shellIndex); //- save the (interpolated) ShellSignalVector as the ith column in the SignalMatrix SignalMatrix.set_column(shellIndex, SignalVector); shellIterator++; shellIndex++; } // ADC averaging Sum(ADC)/n for(unsigned int i = 0 ; i < SignalMatrix.rows(); i++) SignalVector.put(i,SignalMatrix.get_row(i).sum()); - SignalVector/= (double)(m_BValueMap.size()-1); + SignalVector/= (double)(m_B_ValueMap.size()-1); //MITK_INFO << "SignalVector(AVERAGEDADC)"; //MITK_INFO << SignalVector ; // recalculate signal from ADC - calculateSignalFromAdc(SignalVector, m_TargetBValue, BZeroAverage); + calculateSignalFromAdc(SignalVector, m_TargetB_Value, BZeroAverage); //MITK_INFO << "SignalVector(ADC->SIGNAL)"; //MITK_INFO << SignalVector ; for(unsigned int i = 1 ; i < out.Size(); i ++) out.SetElement(i,SignalVector.get(i-1)); //MITK_INFO << "NewSignal: " << out; //MITK_INFO << "-------------"; oit.Set(out); ++oit; ++iit; } } template void MultiShellAdcAverageReconstructionImageFilter ::S_S0Normalization( vnl_vector & vec, const double & S0 ) { for(unsigned int i = 0; i < vec.size(); i++){ vec[i] /= S0; if(vec[i]>1.0) vec[i] = 1; if(vec[i]<0.0) vec[i] = 0; } } template void MultiShellAdcAverageReconstructionImageFilter ::calculateAdcFromSignal( vnl_vector & vec, const double & bValue) { for(unsigned int i = 0; i < vec.size(); i++) vec[i] = log(vec[i])/-bValue; } template void MultiShellAdcAverageReconstructionImageFilter ::calculateSignalFromAdc(vnl_vector & vec, const double & bValue, const double & referenceSignal) { for(unsigned int i = 0 ; i < vec.size(); i++){ //MITK_INFO << vec[i]; //MITK_INFO << bValue; //MITK_INFO << referenceSignal; vec[i] = referenceSignal * exp((-bValue) * vec[i]); } } } // end of namespace diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h index 29f3290768..fc641277ef 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/Reconstruction/itkMultiShellAdcAverageReconstructionImageFilter.h @@ -1,119 +1,119 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _itk_MultiShellAdcAverageReconstructionImageFilter_h_ #define _itk_MultiShellAdcAverageReconstructionImageFilter_h_ #include #include #include namespace itk { /** * \brief Select subset of the input vectors equally distributed over the sphere using an iterative electrostatic repulsion strategy. */ template class MultiShellAdcAverageReconstructionImageFilter : public ImageToImageFilter, itk::VectorImage > { public: typedef MultiShellAdcAverageReconstructionImageFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< itk::VectorImage, itk::VectorImage > Superclass; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(MultiShellAdcAverageReconstructionImageFilter, ImageToImageFilter) typedef TInputScalarType InputScalarType; typedef itk::VectorImage InputImageType; typedef typename InputImageType::PixelType InputPixelType; typedef TOutputScalarType OutputScalarType; typedef itk::VectorImage OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef OutputScalarType BaselineScalarType; typedef BaselineScalarType BaselinePixelType; typedef typename itk::Image BaselineImageType; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef std::vector IndicesVector; typedef std::map BValueMap; GradientDirectionContainerType::Pointer GetOriginalGradientDirections(){return m_OriginalGradientDirections;} void SetOriginalGradientDirections(GradientDirectionContainerType::Pointer ptr){m_OriginalGradientDirections = ptr;} GradientDirectionContainerType::Pointer GetTargetGradientDirections(){return m_TargetGradientDirections;} - double GetTargetBValue(){return m_TargetBValue;} + double GetTargetB_Value(){return m_TargetB_Value;} - double GetBValue(){return m_BValue;} - void SetBValue(double val){m_BValue = val;} + double GetB_Value(){return m_B_Value;} + void SetB_Value(double val){m_B_Value = val;} - void SetOriginalBValueMap(BValueMap inp){m_BValueMap = inp;} + void SetOriginalBValueMap(BValueMap inp){m_B_ValueMap = inp;} protected: MultiShellAdcAverageReconstructionImageFilter(); ~MultiShellAdcAverageReconstructionImageFilter() {} void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType NumberOfThreads ); void S_S0Normalization( vnl_vector & vec, const double & S0 ); void calculateAdcFromSignal( vnl_vector & vec, const double & bValue); void calculateSignalFromAdc( vnl_vector & vec, const double & bValue, const double & referenceSignal); GradientDirectionContainerType::Pointer m_TargetGradientDirections; ///< container for the subsampled output gradient directions GradientDirectionContainerType::Pointer m_OriginalGradientDirections; ///< input gradient directions - BValueMap m_BValueMap; - double m_BValue; + BValueMap m_B_ValueMap; + double m_B_Value; - double m_TargetBValue; + double m_TargetB_Value; std::vector m_WeightsVector; std::vector > m_ShellInterpolationMatrixVector; std::vector m_bZeroIndicesSplitVectors; IndicesVector m_allDirectionsIndicies; unsigned int m_allDirectionsSize; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkMultiShellAdcAverageReconstructionImageFilter.cpp" #endif #endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp index 36b853feb7..50d0df7dcf 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp @@ -1,636 +1,634 @@ /*=================================================================== 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->SetB_Value(m_DiffusionImage->GetB_Value()); filter->Update(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( filter->GetOutput() ); - image->SetB_Value( filter->GetTargetBValue() ); - //image->SetOriginalDirections( filter->GetTargetGradientDirections() ); + image->SetB_Value( filter->GetTargetB_Value() ); 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); m_Controls->m_AdcAverage->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; } }