diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp index 7ca69c9b58..67553574ed 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp @@ -1,796 +1,775 @@ /*========================================================================= 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. =========================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkGibbsTrackingView.h" #include // Qt #include #include #include // MITK #include #include #include #include // ITK #include #include // MISC #include QmitkTrackingWorker::QmitkTrackingWorker(QmitkGibbsTrackingView* view) : m_View(view) { } void QmitkTrackingWorker::run() { m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); MITK_INFO << "Resampling mask images"; // setup resampler typedef itk::ResampleImageFilter ResamplerType; ResamplerType::Pointer resampler = ResamplerType::New(); resampler->SetOutputSpacing( m_View->m_ItkQBallImage->GetSpacing() ); resampler->SetOutputOrigin( m_View->m_ItkQBallImage->GetOrigin() ); resampler->SetOutputDirection( m_View->m_ItkQBallImage->GetDirection() ); resampler->SetSize( m_View->m_ItkQBallImage->GetLargestPossibleRegion().GetSize() ); // resample mask image resampler->SetInput( m_View->m_MaskImage ); resampler->SetDefaultPixelValue(0); resampler->Update(); m_View->m_MaskImage = resampler->GetOutput(); if (m_View->m_GfaImage.IsNotNull()) { ResamplerType::Pointer resampler = ResamplerType::New(); resampler->SetOutputSpacing( m_View->m_ItkQBallImage->GetSpacing() ); resampler->SetOutputOrigin( m_View->m_ItkQBallImage->GetOrigin() ); resampler->SetOutputDirection( m_View->m_ItkQBallImage->GetDirection() ); resampler->SetSize( m_View->m_ItkQBallImage->GetLargestPossibleRegion().GetSize() ); resampler->SetInput( m_View->m_GfaImage ); resampler->SetDefaultPixelValue(0); resampler->Update(); m_View->m_GfaImage = resampler->GetOutput(); } m_View->m_GlobalTracker->SetInput0(m_View->m_ItkQBallImage.GetPointer()); m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); m_View->m_GlobalTracker->SetGfaImage(m_View->m_GfaImage); m_View->m_GlobalTracker->SetTempStart((float)m_View->m_Controls->m_StartTempSlider->value()/100); m_View->m_GlobalTracker->SetTempEnd((float)m_View->m_Controls->m_EndTempSlider->value()/10000); m_View->m_GlobalTracker->SetNumIt(m_View->m_Iterations); m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); m_View->m_GlobalTracker->SetSubtractMean(m_View->m_Controls->m_MeanSubtractionCheckbox->isChecked()); m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); m_View->m_GlobalTracker->SetFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); m_View->m_GlobalTracker->Update(); m_View->m_TrackingThread.quit(); } const std::string QmitkGibbsTrackingView::VIEW_ID = "org.mitk.views.gibbstracking"; QmitkGibbsTrackingView::QmitkGibbsTrackingView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) , m_ThreadIsRunning(false) , m_GlobalTracker(NULL) , m_QBallImage(NULL) , m_MaskImage(NULL) , m_GfaImage(NULL) , m_GfaImageNode(NULL) , m_QBallImageNode(NULL) , m_ItkQBallImage(NULL) , m_FiberBundleNode(NULL) , m_TrackingWorker(this) , m_QBallSelected(false) - , m_FibSelected(false) , m_Iterations(10000000) , m_LastStep(0) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); connect(&m_TrackingThread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } QmitkGibbsTrackingView::~QmitkGibbsTrackingView() { delete m_TrackingTimer; } // update tracking status and generate fiber bundle void QmitkGibbsTrackingView::TimerUpdate() { mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetCurrentStep()-m_LastStep); m_LastStep = m_GlobalTracker->GetCurrentStep(); UpdateTrackingStatus(); GenerateFiberBundle(); } // tell global tractography filter to stop after current step void QmitkGibbsTrackingView::StopGibbsTracking() { if (m_GlobalTracker.IsNull()) return; mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetSteps()-m_LastStep); m_GlobalTracker->SetAbortTracking(true); m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); } // update gui elements and generate fiber bundle after tracking is finished void QmitkGibbsTrackingView::AfterThread() { m_ThreadIsRunning = false; m_TrackingTimer->stop(); UpdateGUI(); UpdateTrackingStatus(); GenerateFiberBundle(); QString paramMessage; if(m_Controls->m_ParticleWeightSlider->value()==0) { m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); paramMessage += "Particle weight was set to " + QString::number(m_GlobalTracker->GetParticleWeight()) + "\n"; } if(m_Controls->m_ParticleWidthSlider->value()==0) { m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); paramMessage += "Particle width was set to " + QString::number(m_GlobalTracker->GetParticleWidth()) + " mm\n"; } if(m_Controls->m_ParticleLengthSlider->value()==0) { m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); paramMessage += "Particle length was set to " + QString::number(m_GlobalTracker->GetParticleLength()) + " mm\n"; } m_FiberBundleNode = NULL; if (paramMessage.length()>0) QMessageBox::information(NULL, "Automatically selected parameters", paramMessage); } // start tracking timer and update gui elements before tracking is started void QmitkGibbsTrackingView::BeforeThread() { m_ThreadIsRunning = true; m_TrackingTime = QTime::currentTime(); m_ElapsedTime = 0; m_TrackingTimer->start(1000); m_LastStep = 0; UpdateGUI(); } // setup gui elements and signal/slot connections void QmitkGibbsTrackingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkGibbsTrackingViewControls; m_Controls->setupUi( parent ); AdvancedSettings(); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); connect( m_Controls->m_SetMaskButton, SIGNAL(clicked()), this, SLOT(SetMask()) ); connect( m_Controls->m_SetGfaButton, SIGNAL(clicked()), this, SLOT(SetGfaImage()) ); connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); connect( m_Controls->m_IterationsSlider, SIGNAL(valueChanged(int)), this, SLOT(SetIterations(int)) ); connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); } } void QmitkGibbsTrackingView::SetInExBalance(int value) { m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); } void QmitkGibbsTrackingView::SetFiberLength(int value) { m_Controls->m_FiberLengthLabel->setText(QString::number(value)); } void QmitkGibbsTrackingView::SetParticleWeight(int value) { if (value>0) { m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); m_Controls->m_GfaFrame->setEnabled(false); } else { m_Controls->m_ParticleWeightLabel->setText("auto"); m_Controls->m_GfaFrame->setEnabled(true); } } void QmitkGibbsTrackingView::SetStartTemp(int value) { m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); } void QmitkGibbsTrackingView::SetEndTemp(int value) { m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); } void QmitkGibbsTrackingView::SetParticleWidth(int value) { if (value>0) m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleWidthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleLength(int value) { if (value>0) m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleLengthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetIterations(int value) { switch(value) { case 0: m_Controls->m_IterationsLabel->setText("Iterations: 10^4"); m_Iterations = 10000; break; case 1: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); m_Iterations = 50000; break; case 2: m_Controls->m_IterationsLabel->setText("Iterations: 10^5"); m_Iterations = 100000; break; case 3: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); m_Iterations = 500000; break; case 4: m_Controls->m_IterationsLabel->setText("Iterations: 10^6"); m_Iterations = 1000000; break; case 5: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); m_Iterations = 5000000; break; case 6: m_Controls->m_IterationsLabel->setText("Iterations: 10^7"); m_Iterations = 10000000; break; case 7: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); m_Iterations = 50000000; break; case 8: m_Controls->m_IterationsLabel->setText("10^8"); m_Iterations = 100000000; break; case 9: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); m_Iterations = 500000000; break; case 10: m_Controls->m_IterationsLabel->setText("Iterations: 10^9"); m_Iterations = 1000000000; break; case 11: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^9"); m_Iterations = 5000000000; break; } } void QmitkGibbsTrackingView::StdMultiWidgetAvailable(QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkGibbsTrackingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } // called if datamanager selection changes void QmitkGibbsTrackingView::OnSelectionChanged( std::vector nodes ) { m_QBallSelected = false; - m_FibSelected = false; // iterate all selected objects for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { m_QBallSelected = true; m_QBallImageNode = node; } - else if (node.IsNotNull() && dynamic_cast(node->GetData())) - { - m_FibSelected = true; - m_FiberBundleNode = node; - } } UpdateGUI(); } // update gui elements displaying trackings status void QmitkGibbsTrackingView::UpdateTrackingStatus() { if (m_GlobalTracker.IsNull()) return; m_ElapsedTime += m_TrackingTime.elapsed()/1000; m_TrackingTime.restart(); unsigned long hours = m_ElapsedTime/3600; unsigned long minutes = (m_ElapsedTime%3600)/60; unsigned long seconds = m_ElapsedTime%60; m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); m_Controls->m_CurrentStepLabel->setText( QString::number(100*(float)m_GlobalTracker->GetCurrentStep()/m_GlobalTracker->GetSteps())+"%" ); m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); } // update gui elements (enable/disable elements and set tooltips) void QmitkGibbsTrackingView::UpdateGUI() { if (!m_ThreadIsRunning && m_QBallSelected) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(true); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_MaskFrame->setEnabled(true); m_Controls->m_GfaFrame->setEnabled(true); m_Controls->m_IterationsSlider->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); m_Controls->m_TrackingStop->setToolTip(""); } else if (!m_ThreadIsRunning) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_MaskFrame->setEnabled(true); m_Controls->m_GfaFrame->setEnabled(true); m_Controls->m_IterationsSlider->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("No Q-Ball image selected."); m_Controls->m_TrackingStop->setToolTip(""); } else { m_Controls->m_TrackingStop->setEnabled(true); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(false); m_Controls->m_MaskFrame->setEnabled(false); m_Controls->m_GfaFrame->setEnabled(false); m_Controls->m_IterationsSlider->setEnabled(false); m_Controls->m_AdvancedFrame->setEnabled(false); m_Controls->m_AdvancedFrame->setVisible(false); m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); } } // show/hide advanced settings frame void QmitkGibbsTrackingView::AdvancedSettings() { m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); } // set mask image data node void QmitkGibbsTrackingView::SetMask() { std::vector nodes = GetDataManagerSelection(); if (nodes.empty()) { m_MaskImageNode = NULL; m_Controls->m_MaskImageEdit->setText("N/A"); return; } for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if (node.IsNotNull() && dynamic_cast(node->GetData())) { m_MaskImageNode = node; m_Controls->m_MaskImageEdit->setText(node->GetName().c_str()); return; } } } // set gfa image data node void QmitkGibbsTrackingView::SetGfaImage() { std::vector nodes = GetDataManagerSelection(); if (nodes.empty()) { m_GfaImageNode = NULL; m_Controls->m_GfaImageEdit->setText("N/A"); return; } for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if (node.IsNotNull() && dynamic_cast(node->GetData())) { m_GfaImageNode = node; m_Controls->m_GfaImageEdit->setText(node->GetName().c_str()); return; } } } // cast image to float template void QmitkGibbsTrackingView::CastToFloat(InputImageType* image, mitk::Image::Pointer outImage) { typedef itk::CastImageFilter ItkCastFilter; typename ItkCastFilter::Pointer itkCaster = ItkCastFilter::New(); itkCaster->SetInput(image); itkCaster->Update(); outImage->InitializeByItk(itkCaster->GetOutput()); outImage->SetVolume(itkCaster->GetOutput()->GetBufferPointer()); } // check for mask and qbi and start tracking thread void QmitkGibbsTrackingView::StartGibbsTracking() { if(m_ThreadIsRunning) { MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; return; } if (!m_QBallSelected) { // Nothing selected. Inform the user and return QMessageBox::information( NULL, "Warning", "Please load and select a qball image before starting image processing."); return; } // a node itself is not very useful, we need its data item (the image) mitk::BaseData* data = m_QBallImageNode->GetData(); if (!data) return; // test if this data item is an image or not (could also be a surface or something totally different) m_QBallImage = dynamic_cast( data ); if (m_QBallImage.IsNull()) return; // cast qbi to itk m_ItkQBallImage = ItkQBallImgType::New(); mitk::CastToItkImage(m_QBallImage, m_ItkQBallImage); // mask image found? // catch exceptions thrown by the itkAccess macros try{ if(m_Controls->m_MaskImageEdit->text().compare("N/A") != 0) { m_MaskImage = 0; if (dynamic_cast(m_MaskImageNode->GetData())) mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); } } catch(...) { QMessageBox::warning(NULL, "Warning", "Incompatible mask image chosen. Processing without masking."); //reset mask image m_MaskImage = NULL; } // if no mask image is selected generate it if( m_MaskImage.IsNull() ) { m_MaskImage = MaskImgType::New(); m_MaskImage->SetSpacing( m_ItkQBallImage->GetSpacing() ); // Set the image spacing m_MaskImage->SetOrigin( m_ItkQBallImage->GetOrigin() ); // Set the image origin m_MaskImage->SetDirection( m_ItkQBallImage->GetDirection() ); // Set the image direction m_MaskImage->SetLargestPossibleRegion( m_ItkQBallImage->GetLargestPossibleRegion()); m_MaskImage->SetBufferedRegion( m_ItkQBallImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); itk::ImageRegionIterator it (m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); for (it = it.Begin(); !it.IsAtEnd(); ++it) { it.Set(1); } } // gfa image found? // catch exceptions thrown by the itkAccess macros try{ if(m_Controls->m_GfaImageEdit->text().compare("N/A") != 0) { m_GfaImage = 0; if (dynamic_cast(m_GfaImageNode->GetData())) mitk::CastToItkImage(dynamic_cast(m_GfaImageNode->GetData()), m_GfaImage); } } catch(...) { QMessageBox::warning(NULL, "Warning", "Incompatible GFA image chosen. Processing without GFA image."); m_GfaImage = NULL; } unsigned int steps = m_Iterations/10000; if (steps<10) steps = 10; m_LastStep = 0; mitk::ProgressBar::GetInstance()->AddStepsToDo(steps); // start worker thread m_TrackingThread.start(QThread::LowestPriority); } // generate mitkFiberBundle from tracking filter output void QmitkGibbsTrackingView::GenerateFiberBundle() { if (m_GlobalTracker.IsNull() || m_ItkQBallImage.IsNull() || m_QBallImage.IsNull() || (!m_Controls->m_VisualizationCheckbox->isChecked() && m_ThreadIsRunning)) return; - m_FiberBundle = mitk::FiberBundle::New(); typedef std::vector< itk::Point > FiberTractType; typedef std::vector< FiberTractType > FiberBundleType; - FiberBundleType* fiberBundle = m_GlobalTracker->GetFiberBundle(); + vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); - for (int i=0; isize(); i++) - { - FiberTractType* tract = &fiberBundle->at(i); - for (int j=0; jsize(); j++) - m_FiberBundle->PushPoint(i, tract->at(j)); - } - m_FiberBundle->initFiberGroup(); - - float bounds[] = {0,0,0}; - bounds[0] = m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(0); - bounds[1] = m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(1); - bounds[2] = m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(2); - m_FiberBundle->SetBounds(bounds); - m_FiberBundle->SetGeometry(m_QBallImage->GetGeometry()); + m_FiberBundle = mitk::FiberBundleX::New(fiberBundle); if (m_FiberBundleNode.IsNotNull()){ GetDefaultDataStorage()->Remove(m_FiberBundleNode); m_FiberBundleNode = 0; } m_FiberBundleNode = mitk::DataNode::New(); m_FiberBundleNode->SetData(m_FiberBundle); QString name(m_QBallImageNode->GetName().c_str()); name += "_FiberBundle"; m_FiberBundleNode->SetName(name.toStdString()); m_FiberBundleNode->SetVisibility(true); if(m_QBallImageNode.IsNull()) GetDataStorage()->Add(m_FiberBundleNode); else GetDataStorage()->Add(m_FiberBundleNode, m_QBallImageNode); } // save current tracking paramters as xml file (.gtp) void QmitkGibbsTrackingView::SaveTrackingParameters() { TiXmlDocument documentXML; TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); documentXML.LinkEndChild( declXML ); TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); mainXML->SetAttribute("file_version", "0.1"); documentXML.LinkEndChild(mainXML); TiXmlElement* paramXML = new TiXmlElement("parameter_set"); paramXML->SetAttribute("iterations", QString::number(m_Iterations).toStdString()); paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); mainXML->LinkEndChild(paramXML); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), QDir::currentPath()+"/param.gtp", tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; if(!filename.endsWith(".gtp")) filename += ".gtp"; documentXML.SaveFile( filename.toStdString() ); } void QmitkGibbsTrackingView::UpdateIteraionsGUI(unsigned long iterations) { switch(iterations) { case 10000: m_Controls->m_IterationsSlider->setValue(0); m_Controls->m_IterationsLabel->setText("Iterations: 10^4"); break; case 50000: m_Controls->m_IterationsSlider->setValue(1); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); break; case 100000: m_Controls->m_IterationsSlider->setValue(2); m_Controls->m_IterationsLabel->setText("Iterations: 10^5"); break; case 500000: m_Controls->m_IterationsSlider->setValue(3); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); break; case 1000000: m_Controls->m_IterationsSlider->setValue(4); m_Controls->m_IterationsLabel->setText("Iterations: 10^6"); break; case 5000000: m_Controls->m_IterationsSlider->setValue(5); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); break; case 10000000: m_Controls->m_IterationsSlider->setValue(6); m_Controls->m_IterationsLabel->setText("Iterations: 10^7"); break; case 50000000: m_Controls->m_IterationsSlider->setValue(7); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); break; case 100000000: m_Controls->m_IterationsSlider->setValue(8); m_Controls->m_IterationsLabel->setText("Iterations: 10^8"); break; case 500000000: m_Controls->m_IterationsSlider->setValue(9); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); break; case 1000000000: m_Controls->m_IterationsSlider->setValue(10); m_Controls->m_IterationsLabel->setText("Iterations: 10^9"); break; case 5000000000: m_Controls->m_IterationsSlider->setValue(11); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^9"); break; } } // load current tracking paramters from xml file (.gtp) void QmitkGibbsTrackingView::LoadTrackingParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; TiXmlDocument doc( filename.toStdString() ); doc.LoadFile(); TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("parameter_set").Element(); QString iterations(pElem->Attribute("iterations")); m_Iterations = iterations.toULong(); UpdateIteraionsGUI(m_Iterations); QString particleLength(pElem->Attribute("particle_length")); float pLength = particleLength.toFloat(); QString particleWidth(pElem->Attribute("particle_width")); float pWidth = particleWidth.toFloat(); if (pLength==0) m_Controls->m_ParticleLengthLabel->setText("auto"); else m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); if (pWidth==0) m_Controls->m_ParticleWidthLabel->setText("auto"); else m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); m_Controls->m_ParticleLengthSlider->setValue(pLength*10); QString partWeight(pElem->Attribute("particle_weight")); m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); m_Controls->m_ParticleWeightLabel->setText(partWeight); QString startTemp(pElem->Attribute("temp_start")); m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); m_Controls->m_StartTempLabel->setText(startTemp); QString endTemp(pElem->Attribute("temp_end")); m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); m_Controls->m_EndTempLabel->setText(endTemp); QString inExBalance(pElem->Attribute("inexbalance")); m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); m_Controls->m_InExBalanceLabel->setText(inExBalance); QString fiberLength(pElem->Attribute("fiber_length")); m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); m_Controls->m_FiberLengthLabel->setText(fiberLength); } diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h index 6f17ec2341..7facb0eed6 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h @@ -1,165 +1,167 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ Version: $Revision: 21975 $ 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. =========================================================================*/ #ifndef QmitkGibbsTrackingView_h #define QmitkGibbsTrackingView_h #include #include #include "ui_QmitkGibbsTrackingViewControls.h" #include #include -#include +#include #include #include +#include +#include class QmitkGibbsTrackingView; class QmitkTrackingWorker : public QObject { Q_OBJECT public: QmitkTrackingWorker(QmitkGibbsTrackingView* view); public slots: void run(); private: QmitkGibbsTrackingView* m_View; }; /*! \brief QmitkGibbsTrackingView \warning This application module is not yet documented. Use "svn blame/praise/annotate" and ask the author to provide basic documentation. \sa QmitkFunctionality \ingroup Functionalities */ typedef itk::Image< float, 3 > FloatImageType; namespace itk { template class GibbsTrackingFilter; } class QmitkGibbsTrackingView : public QmitkFunctionality { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: typedef itk::Image MaskImgType; typedef itk::Vector OdfVectorType; typedef itk::Image ItkQBallImgType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; static const std::string VIEW_ID; QmitkGibbsTrackingView(); virtual ~QmitkGibbsTrackingView(); virtual void CreateQtPartControl(QWidget *parent); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); signals: protected slots: void StartGibbsTracking(); void StopGibbsTracking(); void AfterThread(); void BeforeThread(); void TimerUpdate(); void SetMask(); void SetGfaImage(); void AdvancedSettings(); void SaveTrackingParameters(); void LoadTrackingParameters(); void SetIterations(int value); void SetParticleWidth(int value); void SetParticleLength(int value); void SetInExBalance(int value); void SetFiberLength(int value); void SetParticleWeight(int value); void SetStartTemp(int value); void SetEndTemp(int value); private: // Visualization & GUI void GenerateFiberBundle(); void UpdateGUI(); void UpdateTrackingStatus(); /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); template void CastToFloat(InputImageType* image, typename mitk::Image::Pointer outImage); void UpdateIteraionsGUI(unsigned long iterations); Ui::QmitkGibbsTrackingViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; // data objects - mitk::FiberBundle::Pointer m_FiberBundle; + mitk::FiberBundleX::Pointer m_FiberBundle; MaskImgType::Pointer m_MaskImage; MaskImgType::Pointer m_GfaImage; mitk::QBallImage::Pointer m_QBallImage; ItkQBallImgType::Pointer m_ItkQBallImage; // data nodes mitk::DataNode::Pointer m_QBallImageNode; mitk::DataNode::Pointer m_MaskImageNode; mitk::DataNode::Pointer m_GfaImageNode; mitk::DataNode::Pointer m_FiberBundleNode; // flags etc. bool m_ThreadIsRunning; QTimer* m_TrackingTimer; QTime m_TrackingTime; unsigned long m_ElapsedTime; bool m_QBallSelected; bool m_FibSelected; unsigned long m_Iterations; int m_LastStep; // global tracker and friends itk::SmartPointer m_GlobalTracker; QmitkTrackingWorker m_TrackingWorker; QThread m_TrackingThread; friend class QmitkTrackingWorker; }; #endif // _QMITKGibbsTrackingVIEW_H_INCLUDED diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index 6e9554621a..2e679ea8f2 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,412 +1,471 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ Version: $Revision: 21975 $ 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 "mitkFiberBundleX.h" /* musthave */ //#include // without geometry, fibers are not rendered #include #include #include #include #include #include // baptize array names const char* mitk::FiberBundleX::COLORCODING_ORIENTATION_BASED = "Color_Orient"; const char* mitk::FiberBundleX::COLORCODING_FA_BASED = "Color_FA"; const char* mitk::FiberBundleX::FIBER_ID_ARRAY = "Fiber_IDs"; mitk::FiberBundleX::FiberBundleX(vtkSmartPointer fiberPolyData ) : m_currentColorCoding(NULL) , m_isModified(false) { //generate geometry of passed polydata if (fiberPolyData == NULL) this->m_FiberPolyData = vtkSmartPointer::New(); else this->m_FiberPolyData = fiberPolyData; this->UpdateFiberGeometry(); } mitk::FiberBundleX::~FiberBundleX() { - // Memory Management - m_FiberPolyData->Delete(); + } /* * set computed fibers from tractography algorithms */ void mitk::FiberBundleX::SetFiberPolyData(vtkSmartPointer fiberPD) { if (fiberPD == NULL) this->m_FiberPolyData = vtkSmartPointer::New(); else this->m_FiberPolyData = fiberPD; m_isModified = true; } /* * return fiberbundle as vtkPolyData * Depending on processing of input fibers, this method returns * the latest processed fibers. */ vtkSmartPointer mitk::FiberBundleX::GetFiberPolyData() { return m_FiberPolyData; } /*=================================== *++++ PROCESSING WITH FIBERS +++++++ ====================================*/ void mitk::FiberBundleX::DoColorCodingOrientationbased() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also polydata needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= /* make sure that processing colorcoding is only called when necessary */ if ( m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) && m_FiberPolyData->GetNumberOfPoints() == m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetNumberOfTuples() ) { // fiberstructure is already colorcoded MITK_INFO << " NO NEED TO REGENERATE COLORCODING! " ; return; } /* Finally, execute color calculation */ vtkPoints* extrPoints = m_FiberPolyData->GetPoints(); int numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = sizeof(rgba); vtkUnsignedCharArray * colorsT = vtkUnsignedCharArray::New(); colorsT->Allocate(numOfPoints * componentSize); colorsT->SetNumberOfComponents(componentSize); colorsT->SetName(COLORCODING_ORIENTATION_BASED); /* checkpoint: does polydata contain any fibers */ int numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) { MITK_INFO << "\n ========= Number of Fibers is 0 and below ========= \n"; return; } /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); // MITK_INFO << "Fib#: " << fi << " of " << numOfFibers << " pnts in fiber: " << pointsPerFiber ; /* single fiber checkpoints: is number of points valid */ if (pointsPerFiber > 1) { /* operate on points of single fiber */ for (int i=0; i 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff1[2])); rgba[3] = (unsigned char) (255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff2[2])); rgba[3] = (unsigned char) (255.0); } colorsT->InsertTupleValue(idList[i], rgba); } //end for loop } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; // colorsT->InsertTupleValue(0, rgba); } else { MITK_INFO << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } }//end for loop m_FiberPolyData->GetPointData()->AddArray(colorsT); /*========================= - this is more relevant for renderer than for fiberbundleX datastructure - think about sourcing this to a explicit method which coordinates colorcoding */ this->SetColorCoding(COLORCODING_ORIENTATION_BASED); m_isModified = true; // =========================== //mini test, shall be ported to MITK TESTINGS! if (colorsT->GetSize() != numOfPoints*componentSize) { MITK_INFO << "ALLOCATION ERROR IN INITIATING COLOR ARRAY"; } } void mitk::FiberBundleX::DoGenerateFiberIds() { if (m_FiberPolyData == NULL) return; // for (int i=0; i<10000000; ++i) // { // if(i%500 == 0) // MITK_INFO << i; // } // MITK_INFO << "Generating Fiber Ids"; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInput(m_FiberPolyData); idFiberFilter->CellIdsOn(); // idFiberFilter->PointIdsOn(); // point id's are not needed idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); idFiberFilter->FieldDataOn(); idFiberFilter->Update(); m_FiberIdDataSet = idFiberFilter->GetOutput(); MITK_INFO << "Generating Fiber Ids...[done] | " << m_FiberIdDataSet->GetNumberOfCells(); } void mitk::FiberBundleX::UpdateFiberGeometry() { float min = itk::NumericTraits::min(); float max = itk::NumericTraits::max(); float b[] = {max, min, max, min, max, min}; vtkCellArray* cells = m_FiberPolyData->GetLines(); cells->InitTraversal(); for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p); if (p[0]b[1]) b[1]=p[0]; if (p[1]b[3]) b[3]=p[1]; if (p[2]b[5]) b[5]=p[2]; } } // provide some buffer space at borders for(int i=0; i<=4; i+=2){ b[i] -=10; } for(int i=1; i<=5; i+=2){ b[i] +=10; } mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(true); geometry->SetFloatBounds(b); this->SetGeometry(geometry); } /*============================== *++++ FIBER INFORMATION +++++++ ===============================*/ QStringList mitk::FiberBundleX::GetAvailableColorCodings() { QStringList availableColorCodings; int numColors = m_FiberPolyData->GetPointData()->GetNumberOfArrays(); for(int i=0; iGetPointData()->GetArrayName(i)); } //this controlstructure shall be implemented by the calling method if (availableColorCodings.isEmpty()) MITK_INFO << "no colorcodings available in fiberbundleX"; // for(int i=0; im_currentColorCoding; } void mitk::FiberBundleX::SetColorCoding(const char* requestedColorCoding) { // MITK_INFO << "FbX try to set colorCoding: " << requestedColorCoding << " compare with: " << COLORCODING_ORIENTATION_BASED; if(strcmp (COLORCODING_ORIENTATION_BASED,requestedColorCoding) == 0 ) { this->m_currentColorCoding = (char*) COLORCODING_ORIENTATION_BASED; this->m_isModified = true; } else if(strcmp (COLORCODING_FA_BASED,requestedColorCoding) == 0 ) { this->m_currentColorCoding = (char*) COLORCODING_FA_BASED; this->m_isModified = true; } else { MITK_INFO << "FIBERBUNDLE X: UNKNOWN COLORCODING in FIBERBUNDLEX Datastructure"; this->m_currentColorCoding = "---"; //will cause blank colorcoding of fibers this->m_isModified = true; } } bool mitk::FiberBundleX::isFiberBundleXModified() { return m_isModified; } void mitk::FiberBundleX::setFBXModificationDone() { m_isModified = false; } +void mitk::FiberBundleX::ResampleFibers() +{ + vtkSmartPointer newPoly = vtkSmartPointer::New(); + vtkSmartPointer newCellArray = vtkSmartPointer::New(); + vtkSmartPointer newPoints = vtkSmartPointer::New(); + vtkSmartPointer cellArray = m_FiberPolyData->GetLines(); + for (int i=0; iGetNumberOfLines(); i++) + { + vtkIdType cellId = i; + vtkIdType numPoints = 0; + vtkIdType* points = NULL; + cellArray->GetCell(cellId, numPoints, points); + + vtkSmartPointer container = vtkSmartPointer::New(); + + // insert start + double* point = m_FiberPolyData->GetPoint(points[0]); + vtkIdType pointId = newPoints->InsertNextPoint(point); + container->GetPointIds()->InsertNextId(pointId); + + float Leng = 0; + float dtau = 0; + int cur_p = 1; + int cur_i = 1; + pVector dR; + float normdR; + for (;;) + { + while (dtau <= len && cur_p < numPoints) + { + dR = points[cur_p] - points[cur_p-1]; + normdR = sqrt(dR.norm_square()); + dtau += normdR; + Leng += normdR; + cur_p++; + } + + if (dtau >= len) // if particles reach next voxel + { + point = points[cur_p-1] - dR*( (dtau-len)/normdR ); + pointId = newPoints->InsertNextPoint(point); + container->GetPointIds()->InsertNextId(pointId); + } + else + { + point = m_FiberPolyData->GetPoint(points[numPoints-1]); + pointId = newPoints->InsertNextPoint(point); + container->GetPointIds()->InsertNextId(pointId); + break; + } + dtau = dtau-len; + cur_i++; + } + newCellArray->InsertNextCell(container); + } + newPoly->SetPoints(newPoints); + newPoly->SetLines(newCellArray); + m_FiberPolyData = newPoly; + UpdateFiberGeometry(); +} /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundleX::UpdateOutputInformation() { } void mitk::FiberBundleX::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundleX::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundleX::VerifyRequestedRegion() { return true; } void mitk::FiberBundleX::SetRequestedRegion( itk::DataObject *data ) { } diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h index df312ee528..207d4c3fd3 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h @@ -1,128 +1,129 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Module: $RCSfile$ Language: C++ Date: $Date$ Version: $Revision: 11989 $ 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. =========================================================================*/ /* =============== IMPORTANT TODO =================== * ==== USE vtkSmartPointer<> when necessary ONLY!!!! */ #ifndef _MITK_FiberBundleX_H #define _MITK_FiberBundleX_H //includes for MITK datastructure #include #include "MitkDiffusionImagingExports.h" //includes storing fiberdata #include //may be replaced by class precompile argument #include // may be replaced by class #include // my be replaced by class #include #include namespace mitk { /** * \brief Base Class for Fiber Bundles; */ class MitkDiffusionImaging_EXPORT FiberBundleX : public BaseData { public: - - // names of certain arrays (e.g colorcodings, etc.) + + // names of certain arrays (e.g colorcodings, etc.) static const char* COLORCODING_ORIENTATION_BASED; static const char* COLORCODING_FA_BASED; static const char* FIBER_ID_ARRAY; - + /* friend classes wanna access typedefs ContainerPointType, ContainerTractType, ContainerType */ friend class FiberBundleXWriter; friend class FiberBundleXReader; - + // ======virtual methods must have====== virtual void UpdateOutputInformation(); virtual void SetRequestedRegionToLargestPossibleRegion(); virtual bool RequestedRegionIsOutsideOfTheBufferedRegion(); virtual bool VerifyRequestedRegion(); virtual void SetRequestedRegion( itk::DataObject *data ); //======================================= - + mitkClassMacro( FiberBundleX, BaseData ) itkNewMacro( Self ) //custom constructor with passing argument mitkNewMacro1Param(Self, vtkSmartPointer) - - - + + + /*====FIBERBUNDLE I/O METHODS====*/ void SetFiberPolyData(vtkSmartPointer); //set result of tractography algorithm in vtkPolyData format using vtkPolyLines vtkSmartPointer GetFiberPolyData(); void UpdateFiberGeometry(); char* GetCurrentColorCoding(); QStringList GetAvailableColorCodings(); void SetColorCoding(const char*); - + bool isFiberBundleXModified(); void setFBXModificationDone(); /*===FIBERBUNDLE PROCESSING METHODS====*/ void DoColorCodingOrientationbased(); void DoGenerateFiberIds(); - + void ResampleFibers(); + /*===FIBERBUNDLE ASSESSMENT METHODS====*/ protected: FiberBundleX( vtkSmartPointer fiberPolyData = NULL ); virtual ~FiberBundleX(); private: - + // The following polydata variables are used for fiber- and pointbased representation of the tractography results. As VTK suggests, one vtkPolyData is used to manage vertices and the other for polylines. // FiberPolyData stores all brain fibers using polylines (in world coordinates) // this variable hosts the smoothed fiber data, this data we generate, therefore a smartpointer structure is recommended // vtkSmartPointer m_FiberPolyData; is depricated -// +// // this variable hosts the original fiber data, no smartpointer needed because who or whatever passes this data to FiberBundleX should use vtkSmartPointer structure - + vtkSmartPointer m_FiberPolyData; //this is a common pointer because fiberDataStructure gets passed to this class. m_FiberStructureData is destroyed in the destructor then. - + // this variable contains all additional IDs of Fibers which are needed for efficient fiber manipulation such as extracting etc. vtkSmartPointer m_FiberIdDataSet; char* m_currentColorCoding; //this flag conzerns only visual representation. bool m_isModified; - - + + }; } // namespace mitk #endif /* _MITK_FiberBundleX_H */ diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp index ebe5243717..50b8b03c5f 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp @@ -1,364 +1,155 @@ #ifndef _BUILDFIBRES #define _BUILDFIBRES -//#include "matrix.h" #include #include #include #include using namespace std; #define PI 3.1415926536 - - -#ifdef TIMING - -static struct timeval timeS; - -class PropStats -{ - int N; - int accept; -public: - void clear() { N = 0; accept = 0;} - void propose() {N++;} - void accepted() {accept++;} - - void report(const char *s) - { - mexPrintf("%s #proposals: %8.2fk acceptratio: %.2f \% \n",s,1.0*N/1000.0,100.0*accept/N); - } -}; - - -class Timing -{ -public: - Timing() { time = 0; ncalls = 0;} - void clear() {time = 0; ncalls=0;} - - - long time; - int ncalls; - - void report(const char *s) - { - mexPrintf("%s total: %10.2fms calls: %10.1fk t/call: %10.3fms \n",s,time/1000.0,1.0*ncalls/1000.0,1.0*time/ncalls); - } - - void report_time(const char *s) - { - mexPrintf("%s: %.2fms \n",s,time/1000.0); - } - -}; - -inline void tic(Timing *t) -{ - gettimeofday( &timeS, NULL); - t->time -= (timeS.tv_sec*1000000 + timeS.tv_usec); - t->ncalls++; -} -inline void toc(Timing *t) -{ - gettimeofday( &timeS, NULL); - t->time += (timeS.tv_sec*1000000 + timeS.tv_usec); -} - - -#endif - #include "ParticleGrid.cpp" -class CCAnalysis +#include +#include +#include +#include +#include +#include +#include + +class FiberBuilder { public: Particle *particles; int pcnt; int attrcnt; typedef vector< Particle* > ParticleContainerType; typedef vector< ParticleContainerType* > FiberContainerType; - FiberContainerType* m_FiberContainer; + vtkSmartPointer m_VtkCellArray; + vtkSmartPointer m_VtkPoints; - CCAnalysis(float *points, int numPoints, double spacing[]) + typedef itk::Vector OdfVectorType; + typedef itk::Image ItkQBallImgType; + ItkQBallImgType::Pointer m_ItkQBallImage; + + FiberBuilder(float *points, int numPoints, double spacing[], ItkQBallImgType::Pointer image) { - m_FiberContainer = new FiberContainerType(); + m_ItkQBallImage = image; particles = (Particle*) malloc(sizeof(Particle)*numPoints); pcnt = numPoints; attrcnt = 10; for (int k = 0; k < numPoints; k++) { Particle *p = &(particles[k]); p->R = pVector(points[attrcnt*k]/spacing[0], points[attrcnt*k+1]/spacing[1],points[attrcnt*k+2]/spacing[2]); p->N = pVector(points[attrcnt*k+3],points[attrcnt*k+4],points[attrcnt*k+5]); p->cap = points[attrcnt*k+6]; p->len = points[attrcnt*k+7]; p->mID = (int) points[attrcnt*k+8]; p->pID = (int) points[attrcnt*k+9]; p->ID = k; p->label = 0; } + m_VtkCellArray = vtkSmartPointer::New(); + m_VtkPoints = vtkSmartPointer::New(); } - ~CCAnalysis() + ~FiberBuilder() { - for (int i=0; isize(); i++) - delete(m_FiberContainer->at(i)); - delete(m_FiberContainer); free(particles); } - int iterate(int minSize) + vtkSmartPointer iterate(int minSize) { int cur_label = 1; int numFibers = 0; for (int k = 0; k < pcnt;k++) { Particle *dp = &(particles[k]); if (dp->label == 0) { - ParticleContainerType* container = new ParticleContainerType(); + vtkSmartPointer container = vtkSmartPointer::New(); dp->label = cur_label; dp->numerator = 0; labelPredecessors(dp, container); labelSuccessors(dp, container); - //labelrecursivly(dp, 0); cur_label++; - if(container->size()>minSize){ - m_FiberContainer->push_back(container); + if(container->GetPointIds()->GetNumberOfIds()>minSize) + { + m_VtkCellArray->InsertNextCell(container); numFibers++; } } } - return numFibers; + vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); + fiberPolyData->SetPoints(m_VtkPoints); + fiberPolyData->SetLines(m_VtkCellArray); + return fiberPolyData; + } + + void AddPoint(Particle *dp, vtkSmartPointer container) + { + itk::ContinuousIndex index; + index[0] = dp->R[0]-0.5; + index[1] = dp->R[1]-0.5; + index[2] = dp->R[2]-0.5; + itk::Point point; + m_ItkQBallImage->TransformContinuousIndexToPhysicalPoint( index, point ); + vtkIdType id = m_VtkPoints->InsertNextPoint(point.GetDataPointer()); + container->GetPointIds()->InsertNextId(id); } - void labelPredecessors(Particle *dp, ParticleContainerType* container) + void labelPredecessors(Particle *dp, vtkSmartPointer container) { if (dp->mID != -1 && dp->mID!=dp->ID) { if (dp->ID!=particles[dp->mID].pID) { if (dp->ID==particles[dp->mID].mID) { int tmp = particles[dp->mID].pID; particles[dp->mID].pID = particles[dp->mID].mID; particles[dp->mID].mID = tmp; } } if (particles[dp->mID].label == 0) { particles[dp->mID].label = dp->label; particles[dp->mID].numerator = dp->numerator-1; labelPredecessors(&(particles[dp->mID]), container); } } - container->push_back(dp); + AddPoint(dp, container); } - void labelSuccessors(Particle *dp, ParticleContainerType* container) + + void labelSuccessors(Particle *dp, vtkSmartPointer container) { - if(container->back()->ID!=dp->ID) - container->push_back(dp); + AddPoint(dp, container); if (dp->pID != -1 && dp->pID!=dp->ID) { if (dp->ID!=particles[dp->pID].mID) { if (dp->ID==particles[dp->pID].pID) { int tmp = particles[dp->pID].pID; particles[dp->pID].pID = particles[dp->pID].mID; particles[dp->pID].mID = tmp; } } if (particles[dp->pID].label == 0) { particles[dp->pID].label = dp->label; particles[dp->pID].numerator = dp->numerator+1; labelSuccessors(&(particles[dp->pID]), container); } } } - - void labelrecursivly(Particle *dp, int depth) - { - int label = dp->label; - - if (dp->mID != -1) - { - if (particles[dp->mID].label == 0) - { - particles[dp->mID].label = label; - labelrecursivly(&(particles[dp->mID]),depth+1); - } - } - if (dp->pID != -1) - { - if (particles[dp->pID].label == 0) - { - particles[dp->pID].label = label; - labelrecursivly(&(particles[dp->pID]),depth+1); - } - } - } }; - - - - -//klaus static int cmpfloat2(const void *p1,const void *p2) -//klaus { -//klaus if (((REAL*)p1)[1] > ((REAL*)p2)[1]) -//klaus return 1; -//klaus else -//klaus return -1; -//klaus } - - - -typedef std::vector vecint; - - - -//void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) -//{ -// -// if(nrhs != 2) { -// printf("\nUsage: Df = STderivative(f)\n\n",nrhs); -// printf(" Computes xxx\n"); -// printf(" Parameters:\n"); -// printf(" f - 2D input image of type REAL \n"); -// printf(" Return value Df contains xxx.\n\n"); -// return; -// } else if(nlhs>2) { -// printf("Too many output arguments\n"); -// return; -// } -// -// fprintf(stderr,"building fibers "); fflush(stderr); -// -// int pcnt = 0; -// const mxArray *Points; -// Points = prhs[pcnt++]; -// int numPoints = mxGetN(Points); -// REAL *points = (REAL*) mxGetData(Points); -// -// -// const mxArray *Params; -// Params = prhs[pcnt++]; -// double *params = (double*) mxGetPr(Params); -// -// -// int minnumelements = (int) params[0]; -// -// -// CCAnalysis ccana(points,numPoints); -// -// #ifdef TIMING -// -// #endif -// -// int numc = ccana.iterate(); -// -// fprintf(stderr,"."); fflush(stderr); -// -// vector components(numc); -// -// int i; -// -// for (i = 0;i < ccana.pcnt;i++) -// { -// components[ccana.particles[i].label-1].push_back(i); -// } -// -// fprintf(stderr,"."); fflush(stderr); -// -// for (i = 0; i < numc; i++) -// { -// Particle *last = &(ccana.particles[components[i][0]]); -// last->numerator = 0; -// Particle *next = (last->pID == -1)? 0 : &(ccana.particles[last->pID]); -// for (;;) -// { -// if (next == 0) -// break; -// next->numerator = last->numerator+1; -// int nextID = -1; -// if (next->pID != last->ID) -// nextID = next->pID; -// if (next->mID != last->ID) -// nextID = next->mID; -// last = next; -// next = (nextID == -1)? 0: &(ccana.particles[nextID]); -// if (last->numerator > components[i].size()) // circular -// break; -// } -// -// last = &(ccana.particles[components[i][0]]); -// next = (last->mID == -1)? 0 : &(ccana.particles[last->mID]); -// for (;;) -// { -// if (next == 0) -// break; -// next->numerator = last->numerator-1; -// int nextID = -1; -// if (next->pID != last->ID) -// nextID = next->pID; -// if (next->mID != last->ID) -// nextID = next->mID; -// last = next; -// next = (nextID == -1)? 0: &(ccana.particles[nextID]); -// if (last->numerator < -components[i].size()) // circular -// break; -// } -// -// } -// -// fprintf(stderr,"."); fflush(stderr); -// -// -// #ifdef TIMING -// -// #endif -// -// int index = 0; -// for (i = 0; i < numc; i++) -// { -// if (components[i].size() >= minnumelements) -// { -// index++; -// } -// } -// -// int cdims[] = {index}; -// plhs[0] = mxCreateCellArray(1,cdims); -// -// index = 0; -// for (i = 0; i < numc; i++) -// { -// mxArray *ll = 0; -// if (components[i].size() >= minnumelements) -// { -// ll = mxCreateNumericMatrix(2,components[i].size(),mxGetClassID(Points),mxREAL); -// REAL *dat = (REAL*) mxGetData(ll); -// for (int k = 0; k < components[i].size(); k++) -// { -// dat[2*k] = components[i][k]; -// dat[2*k+1] = ccana.particles[components[i][k]].numerator; -// } -// qsort(dat,components[i].size(),sizeof(REAL)*2,cmpfloat2); -// mxSetCell(plhs[0],index++,ll); -// } -// } -// -// fprintf(stderr,".\n"); fflush(stderr); -// -//} -// - #endif diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index 6ca551d4f2..a0df2f1adc 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,580 +1,546 @@ #include "itkGibbsTrackingFilter.h" #include #include #include "itkPointShell.h" #include "GibbsTracking/BuildFibres.cpp" #pragma GCC visibility push(default) #include #pragma GCC visibility pop #include -#include -#include #include #include -#include #include "GibbsTracking/reparametrize_arclen2.cpp" #include #include #include #include #include #include struct LessDereference { template bool operator()(const T * lhs, const T * rhs) const { return *lhs < *rhs; } }; namespace itk{ template< class TInputOdfImage, class TInputROIImage > GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GibbsTrackingFilter(): m_TempStart(0.1), m_TempEnd(0.001), m_NumIt(500000), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), m_ChempotConnection(10), m_ChempotParticle(0), m_InexBalance(0), m_Chempot2(0.2), m_FiberLength(10), m_AbortTracking(false), m_NumConnections(0), m_NumParticles(0), m_NumAcceptedFibers(0), m_CurrentStep(0), m_SubtractMean(true), m_BuildFibers(false), m_Sampler(NULL), m_Steps(10), m_Memory(0), m_ProposalAcceptance(0), m_GfaImage(NULL) { //this->m_MeasurementFrame.set_identity(); this->SetNumberOfRequiredInputs(2); //Filter needs a DWI image + a Mask Image } template< class TInputOdfImage, class TInputROIImage > GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::~GibbsTrackingFilter(){ delete BESSEL_APPROXCOEFF; if (m_Sampler!=NULL) delete m_Sampler; } template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::ComputeFiberCorrelation(){ // float bD = 15; // vnl_matrix_fixed bDir = // *itk::PointShell >::DistributePointShell(); // const int N = QBALL_ODFSIZE; // vnl_matrix_fixed temp = bDir.transpose(); // vnl_matrix_fixed C = temp*bDir; // vnl_matrix_fixed Q = C; // vnl_vector_fixed mean; // for(int i=0; i repMean; // for (int i=0; i P = Q*Q; // std::vector pointer; // pointer.reserve(N*N); // double * start = C.data_block(); // double * end = start + N*N; // for (double * iter = start; iter != end; ++iter) // { // pointer.push_back(iter); // } // std::sort(pointer.begin(), pointer.end(), LessDereference()); // vnl_vector_fixed alpha; // vnl_vector_fixed beta; // for (int i=0; im_Meanval_sq = (sum*sum)/N; // vnl_vector_fixed alpha_0; // vnl_vector_fixed alpha_2; // vnl_vector_fixed alpha_4; // vnl_vector_fixed alpha_6; // for(int i=0; i T; // T.set_column(0,alpha_0); // T.set_column(1,alpha_2); // T.set_column(2,alpha_4); // T.set_column(3,alpha_6); // vnl_vector_fixed coeff = vnl_matrix_inverse(T).pinverse()*beta; // MITK_INFO << "itkGibbsTrackingFilter: Bessel oefficients: " << coeff; BESSEL_APPROXCOEFF = new float[4]; // BESSEL_APPROXCOEFF[0] = coeff(0); // BESSEL_APPROXCOEFF[1] = coeff(1); // BESSEL_APPROXCOEFF[2] = coeff(2); // BESSEL_APPROXCOEFF[3] = coeff(3); BESSEL_APPROXCOEFF[0] = -0.1714; BESSEL_APPROXCOEFF[1] = 0.5332; BESSEL_APPROXCOEFF[2] = -1.4889; BESSEL_APPROXCOEFF[3] = 2.0389; } // build fibers from tracking result template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::BuildFibers(float* points, int numPoints) { - MITK_INFO << "itkGibbsTrackingFilter: Building fibers ..."; + MITK_INFO << "itkGibbsTrackingFilter: Building fibers ..."; - typename InputQBallImageType::Pointer odfImage - = dynamic_cast(this->GetInput(0)); - double spacing[3]; - spacing[0] = odfImage->GetSpacing().GetElement(0); - spacing[1] = odfImage->GetSpacing().GetElement(1); - spacing[2] = odfImage->GetSpacing().GetElement(2); + double spacing[3]; + spacing[0] = m_ItkQBallImage->GetSpacing().GetElement(0); + spacing[1] = m_ItkQBallImage->GetSpacing().GetElement(1); + spacing[2] = m_ItkQBallImage->GetSpacing().GetElement(2); - // initialize array of particles - CCAnalysis ccana(points, numPoints, spacing); + m_FiberPolyData = FiberPolyDataType::New(); - // label the particles according to fiber affiliation and return number of fibers - int numFibers = ccana.iterate(m_FiberLength); + // initialize array of particles + FiberBuilder fiberBuilder(points, numPoints, spacing, m_ItkQBallImage); + // label the particles according to fiber affiliation and return polydata + m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); + m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); - if (numFibers<=0){ - MITK_INFO << "itkGibbsTrackingFilter: 0 fibers accepted"; - return; - } - - // fill output datastructure - m_FiberBundle.clear(); - for (int i = 0; i < numFibers; i++) - { - vector< Particle* >* particleContainer = ccana.m_FiberContainer->at(i); - - // resample fibers - std::vector< Particle* >* pCon = ResampleFibers(particleContainer, 0.9*spacing[0]); - - FiberTractType tract; - for (int j=0; jsize(); j++) - { - Particle* particle = pCon->at(j); - pVector p = particle->R; - - itk::Point point; - point[0] = p[0]-0.5; - point[1] = p[1]-0.5; - point[2] = p[2]-0.5; - tract.push_back(point); - delete(particle); - } - m_FiberBundle.push_back(tract); - delete(pCon); - } - m_NumAcceptedFibers = numFibers; - MITK_INFO << "itkGibbsTrackingFilter: " << numFibers << " fibers accepted"; + MITK_INFO << "itkGibbsTrackingFilter: " << m_NumAcceptedFibers << " accepted"; } // fill output fiber bundle datastructure template< class TInputOdfImage, class TInputROIImage > -typename GibbsTrackingFilter< TInputOdfImage, TInputROIImage >::FiberBundleType* +typename GibbsTrackingFilter< TInputOdfImage, TInputROIImage >::FiberPolyDataType GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } - return &m_FiberBundle; + return m_FiberPolyData; } // get memory allocated for particle grid template< class TInputOdfImage, class TInputROIImage > float GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GetMemoryUsage() { if (m_Sampler!=NULL) return m_Sampler->m_ParticleGrid.GetMemoryUsage(); return 0; } template< class TInputOdfImage, class TInputROIImage > bool GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::EstimateParticleWeight(){ if (m_GfaImage.IsNull()) { MITK_INFO << "no gfa image found"; return false; } float samplingStart = 1.0; float samplingStop = 0.66; // copy GFA image (original should not be changed) typedef itk::ImageDuplicator< GfaImageType > DuplicateFilterType; DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); duplicator->SetInputImage( m_GfaImage ); duplicator->Update(); m_GfaImage = duplicator->GetOutput(); //// GFA iterator //// typedef ImageRegionIterator< GfaImageType > GfaIteratorType; GfaIteratorType gfaIt(m_GfaImage, m_GfaImage->GetLargestPossibleRegion() ); //// Mask iterator //// typedef ImageRegionConstIterator< MaskImageType > MaskIteratorType; MaskIteratorType maskIt(m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); // set unmasked region of gfa image to 0 gfaIt.GoToBegin(); maskIt.GoToBegin(); while( !gfaIt.IsAtEnd() ) { if(maskIt.Get()<=0) gfaIt.Set(0); ++gfaIt; ++maskIt; } // rescale gfa image to [0,1] typedef itk::RescaleIntensityImageFilter< GfaImageType, GfaImageType > RescaleFilterType; RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New(); rescaleFilter->SetInput( m_GfaImage ); rescaleFilter->SetOutputMaximum( samplingStart ); rescaleFilter->SetOutputMinimum( 0 ); rescaleFilter->Update(); m_GfaImage = rescaleFilter->GetOutput(); gfaIt = GfaIteratorType(m_GfaImage, m_GfaImage->GetLargestPossibleRegion() ); //// Input iterator //// typedef ImageRegionConstIterator< InputQBallImageType > InputIteratorType; InputIteratorType git(m_ItkQBallImage, m_ItkQBallImage->GetLargestPossibleRegion() ); float upper = 0; int count = 0; for(float thr=samplingStart; thr>samplingStop; thr-=0.01) { git.GoToBegin(); gfaIt.GoToBegin(); while( !gfaIt.IsAtEnd() ) { if(gfaIt.Get()>thr) { itk::OrientationDistributionFunction odf(git.Get().GetDataPointer()); upper += odf.GetMaxValue()-odf.GetMeanValue(); ++count; } ++gfaIt; ++git; } } if (count>0) upper /= count; else return false; m_ParticleWeight = upper/6; return true; } // perform global tracking template< class TInputOdfImage, class TInputROIImage > void GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::GenerateData(){ // input qball image m_ItkQBallImage = dynamic_cast(this->GetInput(0)); // approximationscoeffizienten der // teilchenkorrelationen im orientierungsraum // 4er vektor ComputeFiberCorrelation(); // image sizes and spacing int qBallImageSize[4] = {QBALL_ODFSIZE, m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(0), m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(1), m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(2)}; double qBallImageSpacing[3] = {m_ItkQBallImage->GetSpacing().GetElement(0),m_ItkQBallImage->GetSpacing().GetElement(1),m_ItkQBallImage->GetSpacing().GetElement(2)}; // make sure image has enough slices if (qBallImageSize[1]<3 || qBallImageSize[2]<3 || qBallImageSize[3]<3) { MITK_INFO << "itkGibbsTrackingFilter: image size < 3 not supported"; m_AbortTracking = true; } // calculate rotation matrix vnl_matrix_fixed directionMatrix = m_ItkQBallImage->GetDirection().GetVnlMatrix(); vnl_vector_fixed d0 = directionMatrix.get_column(0); d0.normalize(); vnl_vector_fixed d1 = directionMatrix.get_column(1); d1.normalize(); vnl_vector_fixed d2 = directionMatrix.get_column(2); d2.normalize(); directionMatrix.set_column(0, d0); directionMatrix.set_column(1, d1); directionMatrix.set_column(2, d2); vnl_matrix_fixed I = directionMatrix*directionMatrix.transpose(); if(!I.is_identity(mitk::eps)){ MITK_INFO << "itkGibbsTrackingFilter: image direction is not a rotation matrix. Tracking not possible!"; m_AbortTracking = true; } // generate local working copy of image buffer int bufferSize = qBallImageSize[0]*qBallImageSize[1]*qBallImageSize[2]*qBallImageSize[3]; float* qBallImageBuffer = (float*) m_ItkQBallImage->GetBufferPointer(); float* workingQballImage = new float[bufferSize]; for (int i=0; i0 && i%qBallImageSize[0] == 0 && i>0) { sum /= qBallImageSize[0]; for (int j=i-qBallImageSize[0]; jGetBufferPointer(); maskImageSize[0] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(0); maskImageSize[1] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(1); maskImageSize[2] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(2); } else { mask = 0; maskImageSize[0] = qBallImageSize[1]; maskImageSize[1] = qBallImageSize[2]; maskImageSize[2] = qBallImageSize[3]; } int mask_oversamp_mult = maskImageSize[0]/qBallImageSize[1]; // load lookuptable QString applicationDir = QCoreApplication::applicationDirPath(); applicationDir.append("/"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); applicationDir.append("../"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); applicationDir.append("../../"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); std::string lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTBaryCoords.bin"); ifstream BaryCoords; BaryCoords.open(lutPath.c_str(), ios::in | ios::binary); float* coords; if (BaryCoords.is_open()) { float tmp; coords = new float [1630818]; BaryCoords.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { BaryCoords.read((char *)&tmp, sizeof(tmp)); coords[i] = tmp; } BaryCoords.close(); } else { MITK_INFO << "itkGibbsTrackingFilter: unable to open barycoords file"; m_AbortTracking = true; } ifstream Indices; lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTIndices.bin"); Indices.open(lutPath.c_str(), ios::in | ios::binary); int* ind; if (Indices.is_open()) { int tmp; ind = new int [1630818]; Indices.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { Indices.read((char *)&tmp, 4); ind[i] = tmp; } Indices.close(); } else { MITK_INFO << "itkGibbsTrackingFilter: unable to open indices file"; m_AbortTracking = true; } // initialize sphere interpolator with lookuptables SphereInterpolator *sinterp = new SphereInterpolator(coords, ind, QBALL_ODFSIZE, 301, 0.5); // get paramters float minSpacing; if(qBallImageSpacing[0]m_NumIt) { MITK_INFO << "itkGibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } MITK_INFO << "Steps: " << m_Steps; unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); // setup metropolis hastings sampler MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; if (m_Sampler!=NULL) delete m_Sampler; m_Sampler = new RJMCMC(NULL, 0, workingQballImage, qBallImageSize, qBallImageSpacing, cellsize); // setup energy computer MITK_INFO << "itkGibbsTrackingFilter: setting up Energy-computer"; EnergyComputer encomp(workingQballImage,qBallImageSize,qBallImageSpacing,sinterp,&(m_Sampler->m_ParticleGrid),mask,mask_oversamp_mult, directionMatrix); encomp.setParameters(m_ParticleWeight,m_ParticleWidth,m_ChempotConnection*m_ParticleLength*m_ParticleLength,m_ParticleLength,curvatureHardThreshold,m_InexBalance,m_Chempot2); m_Sampler->SetEnergyComputer(&encomp); m_Sampler->SetParameters(m_TempStart,singleIts,m_ParticleLength,curvatureHardThreshold,m_ChempotParticle); // main loop for( int step = 0; step < m_Steps; step++ ) { - if (m_AbortTracking) - break; + if (m_AbortTracking) + break; - m_CurrentStep = step+1; - float temperature = m_TempStart * exp(alpha*(((1.0)*step)/((1.0)*m_Steps))); - MITK_INFO << "itkGibbsTrackingFilter: iterating step " << m_CurrentStep; + m_CurrentStep = step+1; + float temperature = m_TempStart * exp(alpha*(((1.0)*step)/((1.0)*m_Steps))); + MITK_INFO << "itkGibbsTrackingFilter: iterating step " << m_CurrentStep; - m_Sampler->SetTemperature(temperature); - m_Sampler->Iterate(&m_ProposalAcceptance, &m_NumConnections, &m_NumParticles, &m_AbortTracking); + m_Sampler->SetTemperature(temperature); + m_Sampler->Iterate(&m_ProposalAcceptance, &m_NumConnections, &m_NumParticles, &m_AbortTracking); - MITK_INFO << "itkGibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; - MITK_INFO << "itkGibbsTrackingFilter: particles: " << m_NumParticles; - MITK_INFO << "itkGibbsTrackingFilter: connections: " << m_NumConnections; - MITK_INFO << "itkGibbsTrackingFilter: progress: " << 100*(float)step/m_Steps << "%"; + MITK_INFO << "itkGibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; + MITK_INFO << "itkGibbsTrackingFilter: particles: " << m_NumParticles; + MITK_INFO << "itkGibbsTrackingFilter: connections: " << m_NumConnections; + MITK_INFO << "itkGibbsTrackingFilter: progress: " << 100*(float)step/m_Steps << "%"; - if (m_BuildFibers) - { - int numPoints = m_Sampler->m_ParticleGrid.pcnt; - float* points = new float[numPoints*m_Sampler->m_NumAttributes]; - m_Sampler->WriteOutParticles(points); - BuildFibers(points, numPoints); - delete points; - m_BuildFibers = false; - } + if (m_BuildFibers) + { + int numPoints = m_Sampler->m_ParticleGrid.pcnt; + float* points = new float[numPoints*m_Sampler->m_NumAttributes]; + m_Sampler->WriteOutParticles(points); + BuildFibers(points, numPoints); + delete points; + m_BuildFibers = false; + } } int numPoints = m_Sampler->m_ParticleGrid.pcnt; float* points = new float[numPoints*m_Sampler->m_NumAttributes]; m_Sampler->WriteOutParticles(points); BuildFibers(points, numPoints); delete points; delete sinterp; delete coords; delete ind; delete workingQballImage; m_AbortTracking = true; m_BuildFibers = false; MITK_INFO << "itkGibbsTrackingFilter: done generate data"; } } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h index 47400c4542..7c2ecfc200 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h @@ -1,172 +1,176 @@ #ifndef itkGibbsTrackingFilter_h #define itkGibbsTrackingFilter_h #include "itkProcessObject.h" #include "itkVectorContainer.h" #include "itkImage.h" #include "GibbsTracking/pcRJMCMC.cpp" #include "GibbsTracking/auxilary_classes.cpp" #include #include +#include +#include +#include +#include +#include namespace itk{ template< class TInputQBallImage, class TInputROIImage > class GibbsTrackingFilter : public ProcessObject{ public: typedef GibbsTrackingFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; itkNewMacro(Self); itkTypeMacro( GibbsTrackingFilter, ProcessObject ); /** Types for the DWI Input Image **/ typedef TInputQBallImage InputQBallImageType; /** Types for the Mask Image **/ typedef TInputROIImage MaskImageType; typedef typename MaskImageType::Pointer MaskImageTypePointer; - typedef std::vector< itk::Point > FiberTractType; - typedef std::vector< FiberTractType > FiberBundleType; + typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; typedef Image< float, 3 > GfaImageType; typedef typename GfaImageType::Pointer GfaImageTypePointer; itkSetMacro( TempStart, float ); itkGetMacro( TempStart, float ); itkSetMacro( TempEnd, float ); itkGetMacro( TempEnd, float ); itkSetMacro( NumIt, unsigned long ); itkGetMacro( NumIt, unsigned long ); itkSetMacro( ParticleWeight, float ); itkGetMacro( ParticleWeight, float ); /** width of particle sigma (std-dev of gaussian around center) **/ itkSetMacro( ParticleWidth, float ); itkGetMacro( ParticleWidth, float ); /** length of particle from midpoint to ends **/ itkSetMacro( ParticleLength, float ); itkGetMacro( ParticleLength, float ); itkSetMacro( ChempotConnection, float ); itkGetMacro( ChempotConnection, float ); itkSetMacro( ChempotParticle, float ); itkGetMacro( ChempotParticle, float ); itkSetMacro( InexBalance, float ); itkGetMacro( InexBalance, float ); itkSetMacro( Chempot2, float ); itkGetMacro( Chempot2, float ); itkSetMacro( FiberLength, int ); itkGetMacro( FiberLength, int ); itkSetMacro( AbortTracking, bool ); itkGetMacro( AbortTracking, bool ); itkSetMacro( CurrentStep, unsigned long ); itkGetMacro( CurrentStep, unsigned long ); itkSetMacro( SubtractMean, bool); itkGetMacro( SubtractMean, bool); /** Set/Get the Odf Input Image **/ itkSetInputMacro(OdfImage, InputQBallImageType, 0); itkGetInputMacro(OdfImage, InputQBallImageType, 0); /** Set/Get the Input mask image **/ itkSetMacro(MaskImage, MaskImageTypePointer); itkGetMacro(MaskImage, MaskImageTypePointer); itkSetMacro(GfaImage, GfaImageTypePointer); itkGetMacro(GfaImage, GfaImageTypePointer); itkGetMacro(NumParticles, unsigned long); itkGetMacro(NumConnections, unsigned long); itkGetMacro(NumAcceptedFibers, int); itkGetMacro(ProposalAcceptance, float); itkGetMacro(Steps, unsigned int); /** Entry Point For the Algorithm: Is invoked when Update() is called either directly or through itk pipeline propagation **/ void GenerateData(); /** override the Process Object Update because we don't have a dataobject as an outpgnome themeut. We can change this later by wrapping the tractcontainer in a dataobject decorator and letting the Superclass know about it. **/ struct StochasticTractGenerationCallbackStruct{ Pointer Filter; }; virtual void Update(){ this->GenerateData(); } - FiberBundleType* GetFiberBundle(); + FiberPolyDataType GetFiberBundle(); float GetMemoryUsage(); bool EstimateParticleWeight(); protected: GibbsTrackingFilter(); virtual ~GibbsTrackingFilter(); void ComputeFiberCorrelation(); void BuildFibers(float* points, int numPoints); // Input Images typename InputQBallImageType::Pointer m_ItkQBallImage; typename MaskImageType::Pointer m_MaskImage; typename GfaImageType::Pointer m_GfaImage; // Tracking parameters float m_TempStart; // Start temperature float m_TempEnd; // End temperature unsigned long m_NumIt; // Total number of iterations unsigned long m_CurrentStep; // current tracking step float m_ParticleWeight; //w (unitless) float m_ParticleWidth; //sigma (mm) float m_ParticleLength; // ell (mm) float m_ChempotConnection; // gross L (chemisches potential) float m_ChempotParticle;// unbenutzt (immer null, wenn groesser dann insgesamt weniger teilchen) float m_InexBalance; // gewichtung zwischen den lambdas // -5 ... 5 -> nur intern ... nur extern,default 0 float m_Chempot2; // typischerweise 0, // korrektur fuer das geschaetzte integral int m_FiberLength; bool m_AbortTracking; bool m_SubtractMean; int m_NumAcceptedFibers; volatile bool m_BuildFibers; unsigned int m_Steps; float m_Memory; float m_ProposalAcceptance; RJMCMC* m_Sampler; - FiberBundleType m_FiberBundle; + FiberPolyDataType m_FiberPolyData; unsigned long m_NumParticles; unsigned long m_NumConnections; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkGibbsTrackingFilter.cpp" #endif #endif