diff --git a/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.h b/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.h index 0e9cfc8394..be9303c6a4 100644 --- a/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.h +++ b/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.h @@ -1,151 +1,154 @@ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.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_ResidualImageFilter_h_ #define _itk_ResidualImageFilter_h_ #include "itkImageToImageFilter.h" #include namespace itk { template class ResidualImageFilter : public ImageToImageFilter, itk::Image > { public: typedef TInputScalarType InputScalarType; typedef itk::VectorImage InputImageType; typedef typename InputImageType::PixelType InputPixelType; typedef typename InputImageType::RegionType InputImageRegionType; typedef TOutputScalarType OutputScalarType; typedef itk::Image OutputImageType; typedef typename OutputImageType::PixelType OutputPixelType; typedef typename OutputImageType::RegionType OutputImageRegionType; typedef ResidualImageFilter Self; typedef ImageToImageFilter Superclass; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; typedef itk::Image BaselineImageType; itkTypeMacro (ResidualImageFilter, ImageToImageFilter); itkStaticConstMacro (ImageDimension, unsigned int, OutputImageType::ImageDimension); itkNewMacro (Self); void SetSecondDiffusionImage(typename InputImageType::Pointer diffImage) { m_SecondDiffusionImage = diffImage; } std::vector GetQ1() { return m_Q1; } std::vector GetQ3() { return m_Q3; } std::vector GetMeans() { return m_Means; } std::vector GetPercentagesOfOutliers() { return m_PercentagesOfOutliers; } void SetGradients(GradientDirectionContainerType* grads) { m_Gradients = grads; } void SetBaseLineImage(BaselineImageType* baseline) { m_BaseLineImage = baseline; } void SetB0Threshold(InputScalarType threshold) { m_B0Threshold = threshold; } itkSetMacro(B0Index, int) protected: ResidualImageFilter() { m_B0Threshold = 30.0; // default value. allow user to redefine }; ~ResidualImageFilter(){}; void PrintSelf (std::ostream& os, Indent indent) const { Superclass::PrintSelf (os, indent); } void GenerateData(); private: ResidualImageFilter (const Self&); void operator=(const Self&); typename InputImageType::Pointer m_SecondDiffusionImage; std::vector m_Means, m_Q1, m_Q3, m_PercentagesOfOutliers; + // 'Outer' vector: volumes, 'Inner' vector slices + std::vector< std::vector > m_OutliersPerSlice; + GradientDirectionContainerType* m_Gradients; typename BaselineImageType::Pointer m_BaseLineImage; InputScalarType m_B0Threshold; int m_B0Index; }; } // end of namespace #ifndef ITK_MANUAL_INSTANTIATION #include "itkResidualImageFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.txx b/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.txx index 63717b14c9..c351b14f1c 100644 --- a/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.txx +++ b/Modules/DiffusionImaging/Algorithms/itkResidualImageFilter.txx @@ -1,211 +1,217 @@ /*========================================================================= Program: Tensor ToolKit - TTK Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToDiffusionImageFilter.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_ResidualImageFilter_txx_ #define _itk_ResidualImageFilter_txx_ #endif #include "itkResidualImageFilter.h" #include namespace itk { template void ResidualImageFilter ::GenerateData() { typename InputImageType::SizeType size = this->GetInput()->GetLargestPossibleRegion().GetSize(); typename InputImageType::SizeType size2 = m_SecondDiffusionImage->GetLargestPossibleRegion().GetSize(); if(size != size2) { MITK_ERROR << "Sizes do not match"; return; } // Initialize output image typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); outputImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing outputImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin outputImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction outputImage->SetRegions( this->GetInput()->GetLargestPossibleRegion() ); outputImage->Allocate(); outputImage->FillBuffer(0.0); std::vector< std::vector > residuals; + // per slice, per volume + std::vector< std::vector > > residualsPerSlice; // Detrmine number of B0 images int numberB0=0; for(int i=0; iSize(); i++) { GradientDirectionType grad = m_Gradients->ElementAt(i); if(grad[0] < 0.001 && grad[1] < 0.001 && grad[2] < 0.001) { numberB0++; } } residuals.resize(this->GetInput()->GetVectorLength()-numberB0); - - // Calculate the standard residual image and for each volume put all residuals in a vector - for(int x=0; x > sliceResiduals; // residuals per volume for this slice + sliceResiduals.resize(this->GetInput()->GetVectorLength()-numberB0); + for(int y=0; y ix; ix[0] = x; ix[1] = y; ix[2] = z; - typename InputImageType::PixelType p1 = this->GetInput()->GetPixel(ix); typename InputImageType::PixelType p2 = m_SecondDiffusionImage->GetPixel(ix); int s1 = p1.GetSize(); int s2 = p2.GetSize(); if(p1.GetSize() != p2.GetSize()) { MITK_ERROR << "Vector sizes do not match"; return; } if(p1.GetElement(m_B0Index) <= m_B0Threshold) { continue; } double res = 0; int shift = 0; // correction for the skipped B0 images for(int i = 0; iElementAt(i); if(!(grad[0] < 0.001 && grad[1] < 0.001 && grad[2] < 0.001)) { double val1 = (double)p1.GetElement(i); double val2 = (double)p2.GetElement(i); - res += abs(val1-val2); residuals[i-shift].push_back(val1-val2); + sliceResiduals[i-shift].push_back(val1-val2); + } else { shift++; } + } res = res/p1.GetSize(); outputImage->SetPixel(ix, res); } } + + residualsPerSlice.push_back(sliceResiduals); } // for each dw volume: sort the the measured residuals (for each voxel) to enable determining Q1 and Q3; calculate means // determine percentage of errors as described in QUALITY ASSESSMENT THROUGH ANALYSIS OF RESIDUALS OF DIFFUSION IMAGE FITTING // Leemans et al 2008 double q1,q3, median; std::vector< std::vector >::iterator it = residuals.begin(); while(it != residuals.end()) { std::vector res = *it; // sort std::sort(res.begin(), res.end()); q1 = res[0.25*res.size()]; m_Q1.push_back(q1); q3 = res[0.75*res.size()]; m_Q3.push_back(q3); median = res[0.5*res.size()]; double iqr = q3-q1; double outlierThreshold = median + 1.5*iqr; double numberOfOutliers = 0.0; std::vector::iterator resIt = res.begin(); double mean; while(resIt != res.end()) { double f = *resIt; if(f>outlierThreshold) { numberOfOutliers++; } mean += f; ++resIt; } double percOfOutliers = 100 * numberOfOutliers / res.size(); m_PercentagesOfOutliers.push_back(percOfOutliers); mean /= res.size(); m_Means.push_back(mean); ++it; } } } // end of namespace diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp index e088b9809d..238d525c20 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp @@ -1,1076 +1,1128 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date: 2009-05-28 17:19:30 +0200 (Do, 28 Mai 2009) $ Version: $Revision: 17495 $ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html 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. =========================================================================*/ #include "QmitkTensorReconstructionView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include +#include +#include +#include +#include // itk includes #include "itkTimeProbe.h" //#include "itkTensor.h" // mitk includes #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "mitkTeemDiffusionTensor3DReconstructionImageFilter.h" #include "itkDiffusionTensor3DReconstructionImageFilter.h" #include "itkTensorImageToDiffusionImageFilter.h" #include "itkPointShell.h" #include "itkVector.h" #include "itkB0ImageExtractionImageFilter.h" #include "mitkProperties.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include "mitkDiffusionImageMapper.h" #include "mitkLookupTableProperty.h" #include "mitkLookupTable.h" #include "mitkImageStatisticsHolder.h" #include "berryIStructuredSelection.h" #include "berryIWorkbenchWindow.h" #include "berryISelectionService.h" #include #include + const std::string QmitkTensorReconstructionView::VIEW_ID = "org.mitk.views.tensorreconstruction"; #define DI_INFO MITK_INFO("DiffusionImaging") typedef float TTensorPixelType; using namespace berry; struct TrSelListener : ISelectionListener { berryObjectMacro(TrSelListener); TrSelListener(QmitkTensorReconstructionView* view) { m_View = view; } void DoSelectionChanged(ISelection::ConstPointer selection) { // if(!m_View->IsVisible()) // return; // save current selection in member variable m_View->m_CurrentSelection = selection.Cast(); // do something with the selected items if(m_View->m_CurrentSelection) { bool foundDwiVolume = false; bool foundTensorVolume = false; // iterate selection for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin(); i != m_View->m_CurrentSelection->End(); ++i) { // extract datatree node if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); // only look at interesting types if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { foundDwiVolume = true; } // only look at interesting types if(QString("TensorImage").compare(node->GetData()->GetNameOfClass())==0) { foundTensorVolume = true; } } } m_View->m_Controls->m_ItkReconstruction->setEnabled(foundDwiVolume); m_View->m_Controls->m_TeemReconstruction->setEnabled(foundDwiVolume); m_View->m_Controls->m_TensorsToDWIButton->setEnabled(foundTensorVolume); m_View->m_Controls->m_TensorsToQbiButton->setEnabled(foundTensorVolume); m_View->m_Controls->m_ResidualButton->setEnabled(foundDwiVolume && foundTensorVolume); m_View->m_Controls->m_PercentagesOfOutliers->setEnabled(foundDwiVolume && foundTensorVolume); } } void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection) { // check, if selection comes from datamanager if (part) { QString partname(part->GetPartName().c_str()); if(partname.compare("Datamanager")==0) { // apply selection DoSelectionChanged(selection); } } } QmitkTensorReconstructionView* m_View; }; QmitkTensorReconstructionView::QmitkTensorReconstructionView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL) { + + + + } QmitkTensorReconstructionView::QmitkTensorReconstructionView(const QmitkTensorReconstructionView& other) { Q_UNUSED(other) throw std::runtime_error("Copy constructor not implemented"); } QmitkTensorReconstructionView::~QmitkTensorReconstructionView() { this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->RemovePostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener); } void QmitkTensorReconstructionView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkTensorReconstructionViewControls; m_Controls->setupUi(parent); this->CreateConnections(); QStringList items; items << "LLS (Linear Least Squares)" << "MLE (Maximum Likelihood)" << "NLS (Nonlinear Least Squares)" << "WLS (Weighted Least Squares)"; m_Controls->m_TensorEstimationTeemEstimationMethodCombo->addItems(items); m_Controls->m_TensorEstimationTeemEstimationMethodCombo->setCurrentIndex(0); m_Controls->m_TensorEstimationManualThreashold->setChecked(false); m_Controls->m_TensorEstimationTeemSigmaEdit->setText("NaN"); m_Controls->m_TensorEstimationTeemNumItsSpin->setValue(1); m_Controls->m_TensorEstimationTeemFuzzyEdit->setText("0.0"); m_Controls->m_TensorEstimationTeemMinValEdit->setText("1.0"); m_Controls->m_TensorEstimationTeemNumItsLabel_2->setEnabled(true); m_Controls->m_TensorEstimationTeemNumItsSpin->setEnabled(true); m_Controls->m_TensorsToDWIBValueEdit->setText("1000"); Advanced1CheckboxClicked(); Advanced2CheckboxClicked(); TeemCheckboxClicked(); #ifndef DIFFUSION_IMAGING_EXTENDED m_Controls->m_TeemToggle->setVisible(false); #endif // define data type for combobox //m_Controls->m_ImageSelector->SetDataStorage( this->GetDefaultDataStorage() ); //m_Controls->m_ImageSelector->SetPredicate( mitk::NodePredicateDataType::New("DiffusionImage") ); } m_SelListener = berry::ISelectionListener::Pointer(new TrSelListener(this)); this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->AddPostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener); berry::ISelection::ConstPointer sel( this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); m_CurrentSelection = sel.Cast(); m_SelListener.Cast()->DoSelectionChanged(sel); + + + + + + // Create some QImage + QImage image(3, 3, QImage::Format_Indexed8); + QRgb value; + + value = qRgb(122, 163, 39); // 0xff7aa327 + image.setColor(0, value); + + value = qRgb(237, 187, 51); // 0xffedba31 + image.setColor(1, value); + value = qRgb(189, 149, 39); // 0xffbd9527 + + image.setColor(2, value); + + image.setPixel(0, 1, 0); + image.setPixel(1, 0, 0); + image.setPixel(1, 1, 2); + image.setPixel(2, 1, 1); + + int dotX = image.dotsPerMeterX(); + int dotY= image.dotsPerMeterY(); + + image.setDotsPerMeterX(10*dotX); + image.setDotsPerMeterY(10*dotY); + + QGraphicsScene* scene = new QGraphicsScene; + + QGraphicsPixmapItem *item = new QGraphicsPixmapItem( QPixmap::fromImage(image), 0, scene); + + item->setScale(25); + + + //QGraphicsView view( &scene ); + m_Controls->m_PerSliceView->setRenderHints( QPainter::Antialiasing ); + + m_Controls->m_PerSliceView->setScene(scene); + m_Controls->m_PerSliceView->show(); + + m_Controls->m_PerSliceView->repaint(); } void QmitkTensorReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { berry::ISelection::ConstPointer sel( this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); m_CurrentSelection = sel.Cast(); m_SelListener.Cast()->DoSelectionChanged(sel); m_MultiWidget = &stdMultiWidget; } void QmitkTensorReconstructionView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkTensorReconstructionView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_TeemToggle), SIGNAL(clicked()), this, SLOT(TeemCheckboxClicked()) ); connect( (QObject*)(m_Controls->m_ItkReconstruction), SIGNAL(clicked()), this, SLOT(ItkReconstruction()) ); connect( (QObject*)(m_Controls->m_TeemReconstruction), SIGNAL(clicked()), this, SLOT(TeemReconstruction()) ); connect( (QObject*)(m_Controls->m_TensorEstimationTeemEstimationMethodCombo), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) ); connect( (QObject*)(m_Controls->m_Advanced1), SIGNAL(clicked()), this, SLOT(Advanced1CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_Advanced2), SIGNAL(clicked()), this, SLOT(Advanced2CheckboxClicked()) ); connect( (QObject*)(m_Controls->m_TensorEstimationManualThreashold), SIGNAL(clicked()), this, SLOT(ManualThresholdClicked()) ); connect( (QObject*)(m_Controls->m_TensorsToDWIButton), SIGNAL(clicked()), this, SLOT(TensorsToDWI()) ); connect( (QObject*)(m_Controls->m_TensorsToQbiButton), SIGNAL(clicked()), this, SLOT(TensorsToQbi()) ); connect( (QObject*)(m_Controls->m_ResidualButton), SIGNAL(clicked()), this, SLOT(ResidualCalculation()) ); } } void QmitkTensorReconstructionView::TeemCheckboxClicked() { m_Controls->groupBox_3->setVisible(m_Controls-> m_TeemToggle->isChecked()); } void QmitkTensorReconstructionView::Advanced1CheckboxClicked() { bool check = m_Controls-> m_Advanced1->isChecked(); m_Controls->frame->setVisible(check); } void QmitkTensorReconstructionView::Advanced2CheckboxClicked() { bool check = m_Controls-> m_Advanced2->isChecked(); m_Controls->frame_2->setVisible(check); } void QmitkTensorReconstructionView::ManualThresholdClicked() { m_Controls->m_TensorReconstructionThreasholdEdit_2->setEnabled( m_Controls->m_TensorEstimationManualThreashold->isChecked()); } void QmitkTensorReconstructionView::Activated() { QmitkFunctionality::Activated(); } void QmitkTensorReconstructionView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkTensorReconstructionView::MethodChoosen(int method) { m_Controls->m_TensorEstimationTeemNumItsLabel_2->setEnabled(method==3); m_Controls->m_TensorEstimationTeemNumItsSpin->setEnabled(method==3); } void QmitkTensorReconstructionView::ResidualCalculation() { // Extract dwi and dti from current selection // In case of multiple selections, take the first one, since taking all combinations is not meaningful if(m_CurrentSelection) { mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); mitk::DiffusionImage::Pointer diffImage = mitk::DiffusionImage::New(); typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer tensorImage; std::string nodename; for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); i != m_CurrentSelection->End(); ++i) { if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { diffImage = static_cast*>((node)->GetData()); } else if((QString("TensorImage").compare(node->GetData()->GetNameOfClass())==0)) { mitk::TensorImage* mitkVol; mitkVol = static_cast((node)->GetData()); mitk::CastToItkImage(mitkVol, tensorImage); node->GetStringProperty("name", nodename); } } } typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; FilterType::GradientListType gradientList; mitk::DiffusionImage::GradientDirectionContainerType* gradients = diffImage->GetDirections(); // Copy gradients vectors from gradients to gradientList for(int i=0; iSize(); i++) { mitk::DiffusionImage::GradientDirectionType vec = gradients->at(i); itk::Vector grad; grad[0] = vec[0]; grad[1] = vec[1]; grad[2] = vec[2]; gradientList.push_back(grad); } // Find the min and the max values from a baseline image mitk::ImageStatisticsHolder *stats = diffImage->GetStatistics(); //Initialize filter that calculates the modeled diffusion weighted signals FilterType::Pointer filter = FilterType::New(); filter->SetInput( tensorImage ); filter->SetBValue(diffImage->GetB_Value()); filter->SetGradientList(gradientList); filter->SetMin(stats->GetScalarValueMin()); filter->SetMax(500); filter->Update(); // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(diffImage->GetB_Value()); image->SetDirections(gradientList); image->SetOriginalDirections(gradientList); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dwi"); node->SetName(newname.toAscii()); GetDefaultDataStorage()->Add(node); std::vector b0Indices = image->GetB0Indices(); /* // Extract B0 typedef itk::B0ImageExtractionImageFilter BaslineFilterType; BaslineFilterType::Pointer baseFilter = BaslineFilterType::New(); baseFilter->SetInput(image->GetVectorImage()); baseFilter->SetDirections(image->GetDirections()); baseFilter->Update(); mitk::Image::Pointer boImage = mitk::Image::New(); boImage->InitializeByItk( baseFilter->GetOutput() ); boImage->SetVolume( baseFilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer b0Node=mitk::DataNode::New(); b0Node->SetData( boImage ); b0Node->SetProperty( "name", mitk::StringProperty::New("baseline")); GetDefaultDataStorage()->Add(b0Node); */ typedef itk::ResidualImageFilter ResidualImageFilterType; ResidualImageFilterType::Pointer residualFilter = ResidualImageFilterType::New(); residualFilter->SetInput(diffImage->GetVectorImage()); residualFilter->SetSecondDiffusionImage(image->GetVectorImage()); residualFilter->SetGradients(gradients); residualFilter->SetB0Index(b0Indices[0]); residualFilter->SetB0Threshold(30); residualFilter->Update(); itk::Image::Pointer residualImage = itk::Image::New(); residualImage = residualFilter->GetOutput(); mitk::Image::Pointer mitkResImg = mitk::Image::New(); mitk::CastToMitkImage(residualImage, mitkResImg); stats = mitkResImg->GetStatistics(); float min = stats->GetScalarValueMin(); float max = stats->GetScalarValueMax(); mitk::LookupTableProperty::Pointer lutProp = mitk::LookupTableProperty::New(); mitk::LookupTable::Pointer lut = mitk::LookupTable::New(); vtkSmartPointer lookupTable = vtkSmartPointer::New(); lookupTable->SetTableRange(min, max); // If you don't want to use the whole color range, you can use // SetValueRange, SetHueRange, and SetSaturationRange lookupTable->Build(); int size = lookupTable->GetTable()->GetSize(); vtkSmartPointer reversedlookupTable = vtkSmartPointer::New(); reversedlookupTable->SetTableRange(min+1, max); reversedlookupTable->Build(); for(int i=0; i<256; i++) { double* rgba = reversedlookupTable->GetTableValue(255-i); lookupTable->SetTableValue(i, rgba[0], rgba[1], rgba[2], rgba[3]); } lut->SetVtkLookupTable(lookupTable); lutProp->SetLookupTable(lut); // Create lookuptable mitk::DataNode::Pointer resNode=mitk::DataNode::New(); resNode->SetData( mitkResImg ); resNode->SetName("Residual Image"); resNode->SetProperty("LookupTable", lutProp); bool b; resNode->GetBoolProperty("use color", b); resNode->SetBoolProperty("use color", false); GetDefaultDataStorage()->Add(resNode); m_MultiWidget->RequestUpdate(); // Draw Graph std::vector means = residualFilter->GetMeans(); std::vector q1s = residualFilter->GetQ1(); std::vector q3s = residualFilter->GetQ3(); std::vector percentagesOfOUtliers = residualFilter->GetPercentagesOfOutliers(); m_Controls->m_ResidualAnalysis->SetMeans(means); m_Controls->m_ResidualAnalysis->SetQ1(q1s); m_Controls->m_ResidualAnalysis->SetQ3(q3s); m_Controls->m_ResidualAnalysis->SetPercentagesOfOutliers(percentagesOfOUtliers); if(m_Controls->m_PercentagesOfOutliers->isChecked()) { m_Controls->m_ResidualAnalysis->DrawPercentagesOfOutliers(); } else { m_Controls->m_ResidualAnalysis->DrawMeans(); } } } void QmitkTensorReconstructionView::ItkReconstruction() { Reconstruct(0); } void QmitkTensorReconstructionView::TeemReconstruction() { Reconstruct(1); } void QmitkTensorReconstructionView::Reconstruct(int method) { if (m_CurrentSelection) { mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); int at = 0; for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); i != m_CurrentSelection->End(); ++i) { if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0) { set->InsertElement(at++, node); } } } if(method == 0) { ItkTensorReconstruction(set); } if(method == 1) { TeemTensorReconstruction(set); } } } void QmitkTensorReconstructionView::ItkTensorReconstruction (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // TENSOR RECONSTRUCTION clock.Start(); MBI_INFO << "Tensor reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef itk::DiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); tensorReconstructionFilter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() ); tensorReconstructionFilter->SetBValue(vols->GetB_Value()); tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreasholdEdit->text().toFloat() ); tensorReconstructionFilter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s."; // TENSORS TO DATATREE mitk::TensorImage::Pointer image = mitk::TensorImage::New(); typedef itk::Image, 3> TensorImageType; TensorImageType::Pointer tensorImage; tensorImage = tensorReconstructionFilter->GetOutput(); // Check the tensor for negative eigenvalues if(m_Controls->m_CheckNegativeEigenvalues->isChecked()) { typedef itk::ImageRegionIterator TensorImageIteratorType; TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetRequestedRegion()); tensorIt.GoToBegin(); while(!tensorIt.IsAtEnd()) { typedef itk::DiffusionTensor3D TensorType; //typedef itk::Tensor TensorType2; TensorType tensor = tensorIt.Get(); // TensorType2 tensor2; /* for(int i=0; i SymEigenSystemType; SymEigenSystemType eig (tensor2.GetVnlMatrix()); for(unsigned int i=0; iInitializeByItk( tensorImage.GetPointer() ); image->SetVolume( tensorReconstructionFilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dti"); SetDefaultNodeProperties(node, newname.toStdString()); nodes.push_back(node); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } void QmitkTensorReconstructionView::TeemTensorReconstruction (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); ++itemiter; // TENSOR RECONSTRUCTION clock.Start(); MBI_INFO << "Teem Tensor reconstruction "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "Teem Tensor reconstruction for %s", nodename.c_str()).toAscii()); typedef mitk::TeemDiffusionTensor3DReconstructionImageFilter< DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType; TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter = TensorReconstructionImageFilterType::New(); tensorReconstructionFilter->SetInput( vols ); if(!m_Controls->m_TensorEstimationTeemSigmaEdit->text().contains(QString("NaN"))) tensorReconstructionFilter->SetSigma( m_Controls->m_TensorEstimationTeemSigmaEdit->text().toFloat() ); switch(m_Controls->m_TensorEstimationTeemEstimationMethodCombo->currentIndex()) { // items << "LLS (Linear Least Squares)" //<< "MLE (Maximum Likelihood)" //<< "NLS (Nonlinear Least Squares)" //<< "WLS (Weighted Least Squares)"; case 0: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsLLS); break; case 1: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsMLE); break; case 2: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsNLS); break; case 3: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsWLS); break; default: tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsLLS); } tensorReconstructionFilter->SetNumIterations( m_Controls->m_TensorEstimationTeemNumItsSpin->value() ); if(m_Controls->m_TensorEstimationManualThreashold->isChecked()) tensorReconstructionFilter->SetConfidenceThreshold( m_Controls->m_TensorReconstructionThreasholdEdit_2->text().toDouble() ); tensorReconstructionFilter->SetConfidenceFuzzyness( m_Controls->m_TensorEstimationTeemFuzzyEdit->text().toFloat() ); tensorReconstructionFilter->SetMinPlausibleValue( m_Controls->m_TensorEstimationTeemMinValEdit->text().toDouble() ); tensorReconstructionFilter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ; // TENSORS TO DATATREE mitk::DataNode::Pointer node2=mitk::DataNode::New(); node2->SetData( tensorReconstructionFilter->GetOutputItk() ); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dtix"); SetDefaultNodeProperties(node2, newname.toStdString()); nodes.push_back(node2); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } void QmitkTensorReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name) { node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) ); node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) ); node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New()); node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New()); node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2)); node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1)); node->SetProperty( "visible", mitk::BoolProperty::New( true ) ); node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) ); node->SetProperty ("layer", mitk::IntProperty::New(100)); node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) ); //node->SetProperty( "opacity", mitk::FloatProperty::New(1.0f) ); node->SetProperty( "name", mitk::StringProperty::New(name) ); } //node->SetProperty( "volumerendering", mitk::BoolProperty::New( false ) ); //node->SetProperty( "use color", mitk::BoolProperty::New( true ) ); //node->SetProperty( "texture interpolation", mitk::BoolProperty::New( true ) ); //node->SetProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() ); //node->SetProperty( "layer", mitk::IntProperty::New(0)); //node->SetProperty( "in plane resample extent by geometry", mitk::BoolProperty::New( false ) ); //node->SetOpacity(1.0f); //node->SetColor(1.0,1.0,1.0); //node->SetVisibility(true); //node->SetProperty( "IsTensorVolume", mitk::BoolProperty::New( true ) ); //mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); //mitk::LevelWindow levelwindow; //// levelwindow.SetAuto( image ); //levWinProp->SetLevelWindow( levelwindow ); //node->GetPropertyList()->SetProperty( "levelwindow", levWinProp ); //// add a default rainbow lookup table for color mapping //if(!node->GetProperty("LookupTable")) //{ // mitk::LookupTable::Pointer mitkLut = mitk::LookupTable::New(); // vtkLookupTable* vtkLut = mitkLut->GetVtkLookupTable(); // vtkLut->SetHueRange(0.6667, 0.0); // vtkLut->SetTableRange(0.0, 20.0); // vtkLut->Build(); // mitk::LookupTableProperty::Pointer mitkLutProp = mitk::LookupTableProperty::New(); // mitkLutProp->SetLookupTable(mitkLut); // node->SetProperty( "LookupTable", mitkLutProp ); //} //if(!node->GetProperty("binary")) // node->SetProperty( "binary", mitk::BoolProperty::New( false ) ); //// add a default transfer function //mitk::TransferFunction::Pointer tf = mitk::TransferFunction::New(); //node->SetProperty ( "TransferFunction", mitk::TransferFunctionProperty::New ( tf.GetPointer() ) ); //// set foldername as string property //mitk::StringProperty::Pointer nameProp = mitk::StringProperty::New( name ); //node->SetProperty( "name", nameProp ); void QmitkTensorReconstructionView::TensorsToDWI() { if (m_CurrentSelection) { mitk::DataStorage::SetOfObjects::Pointer set = mitk::DataStorage::SetOfObjects::New(); int at = 0; for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); i != m_CurrentSelection->End(); ++i) { if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) { mitk::DataNode::Pointer node = nodeObj->GetDataNode(); if(QString("TensorImage").compare(node->GetData()->GetNameOfClass())==0) { set->InsertElement(at++, node); } } } DoTensorsToDWI(set); } } void QmitkTensorReconstructionView::TensorsToQbi() { std::vector nodes = this->GetDataManagerSelection(); for (int i=0; i TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(dynamic_cast(tensorImageNode->GetData()), itkvol); typedef itk::TensorImageToQBallImageFilter< TTensorPixelType, TTensorPixelType > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->Update(); typedef itk::Vector OutputPixelType; typedef itk::Image OutputImageType; mitk::QBallImage::Pointer image = mitk::QBallImage::New(); OutputImageType::Pointer outimg = filter->GetOutput(); image->InitializeByItk( outimg.GetPointer() ); image->SetVolume( outimg->GetBufferPointer() ); mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData( image ); QString newname; newname = newname.append(tensorImageNode->GetName().c_str()); newname = newname.append("_qbi"); node->SetName(newname.toAscii()); GetDefaultDataStorage()->Add(node); } } void QmitkTensorReconstructionView::OnSelectionChanged( std::vector nodes ) { if ( !this->IsVisible() ) return; } template std::vector > QmitkTensorReconstructionView::MakeGradientList() { std::vector > retval; vnl_matrix_fixed* U = itk::PointShell >::DistributePointShell(); for(int i=0; i v; v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i); retval.push_back(v); } // Add 0 vector for B0 itk::Vector v; v.Fill(0.0); retval.push_back(v); return retval; } void QmitkTensorReconstructionView::DoTensorsToDWI (mitk::DataStorage::SetOfObjects::Pointer inImages) { try { itk::TimeProbe clock; int nrFiles = inImages->size(); if (!nrFiles) return; QString status; mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles); mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { std::string nodename; (*itemiter)->GetStringProperty("name", nodename); mitk::TensorImage* vol = static_cast((*itemiter)->GetData()); ++itemiter; typedef float TTensorPixelType; typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType; typedef itk::Image< TensorPixelType, 3 > TensorImageType; TensorImageType::Pointer itkvol = TensorImageType::New(); mitk::CastToItkImage(vol, itkvol); typedef itk::TensorImageToDiffusionImageFilter< TTensorPixelType, DiffusionPixelType > FilterType; FilterType::GradientListType gradientList; switch(m_Controls->m_TensorsToDWINumDirsSelect->currentIndex()) { case 0: gradientList = MakeGradientList<12>(); break; case 1: gradientList = MakeGradientList<42>(); break; case 2: gradientList = MakeGradientList<92>(); break; case 3: gradientList = MakeGradientList<162>(); break; case 4: gradientList = MakeGradientList<252>(); break; case 5: gradientList = MakeGradientList<362>(); break; case 6: gradientList = MakeGradientList<492>(); break; case 7: gradientList = MakeGradientList<642>(); break; case 8: gradientList = MakeGradientList<812>(); break; case 9: gradientList = MakeGradientList<1002>(); break; default: gradientList = MakeGradientList<92>(); } double bVal = m_Controls->m_TensorsToDWIBValueEdit->text().toDouble(); // DWI ESTIMATION clock.Start(); MBI_INFO << "DWI Estimation "; mitk::StatusBar::GetInstance()->DisplayText(status.sprintf( "DWI Estimation for %s", nodename.c_str()).toAscii()); FilterType::Pointer filter = FilterType::New(); filter->SetInput( itkvol ); filter->SetBValue(bVal); filter->SetGradientList(gradientList); //filter->SetNumberOfThreads(1); filter->Update(); clock.Stop(); MBI_DEBUG << "took " << clock.GetMeanTime() << "s."; itk::Vector v; v[0] = 0; v[1] = 0; v[2] = 0; gradientList.push_back(v); // TENSORS TO DATATREE mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(bVal); image->SetDirections(gradientList); image->SetOriginalDirections(gradientList); image->InitializeFromVectorImage(); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( image ); mitk::DiffusionImageMapper::SetDefaultProperties(node); QString newname; newname = newname.append(nodename.c_str()); newname = newname.append("_dwi"); node->SetName(newname.toAscii()); nodes.push_back(node); mitk::ProgressBar::GetInstance()->Progress(); } std::vector::iterator nodeIt; for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt) GetDefaultDataStorage()->Add(*nodeIt); mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii()); m_MultiWidget->RequestUpdate(); } catch (itk::ExceptionObject &ex) { MBI_INFO << ex ; return ; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui index 26a380ab64..fac9d3fa3e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui @@ -1,598 +1,622 @@ QmitkTensorReconstructionViewControls 0 0 345 1303 0 0 true QmitkTensorReconstructionViewControls ITK Reconstruction Advanced Settings false QFrame::StyledPanel QFrame::Raised QFormLayout::AllNonFixedFieldsGrow B0 Threshold false 0 Check for negative eigenvalues QFrame::StyledPanel QFrame::Raised false ITK Tensor Reconstruction Estimate Diffusion Image from Tensors QFrame::StyledPanel QFrame::Raised QFormLayout::AllNonFixedFieldsGrow 6 6 9 how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" B-Value false 0 0 how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" # Gradient Directions 3 12 42 92 162 252 362 492 642 812 1002 false Diffusion Image Estimation Estimate Q-Ball Image from Tensors false Calculate ODF value as tensor value in the according direction Q-Ball Image Estimation true Teem Reconstruction true Teem Reconstruction Advanced Settings QFrame::StyledPanel QFrame::Raised QFormLayout::AllNonFixedFieldsGrow important in case of method wls # Iterations false how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" Fuzzy confidence false 0 0 how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0" minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0" minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0" minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0" Min plausible value false Rician noise parameter (float) Rician noise parameter (float) Rician noise parameter (float) Sigma false 0 0 minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0" minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0" minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0" 0 0 Rician noise parameter (float) Rician noise parameter (float) Rician noise parameter (float) Method B0-Threshold false 0 false Teem Tensor Reconstruction Residuals true true QFrame::StyledPanel QFrame::Raised false Calculate the residual from a dti and a dwi iimage Residual Image Calculation - - - - 200 - 300 - - + + + 1 + + + + Per volume + + + + + + + 200 + 300 + + + + + + + + + Per slice + + + + + + + false percentages of error Qt::Vertical 20 1150 QmitkResidualAnalysisWidget QWidget
QmitkResidualAnalysisWidget.h
1