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..25914ced2c 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,786 @@ /*========================================================================= 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); + + double qBallImageSpacing[3] = {m_ItkQBallImage->GetSpacing().GetElement(0),m_ItkQBallImage->GetSpacing().GetElement(1),m_ItkQBallImage->GetSpacing().GetElement(2)}; + float minSpacing; + if(qBallImageSpacing[0]ResampleFibers(minSpacing); 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..837669f448 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,412 +1,487 @@ /*========================================================================= 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; } +// Resample fiber to get equidistant points +void mitk::FiberBundleX::ResampleFibers(float len) +{ + vtkSmartPointer newPoly = vtkSmartPointer::New(); + vtkSmartPointer newCellArray = vtkSmartPointer::New(); + vtkSmartPointer newPoints = vtkSmartPointer::New(); + + vtkSmartPointer vLines = m_FiberPolyData->GetLines(); + vLines->InitTraversal(); + int numberOfLines = vLines->GetNumberOfCells(); + + for (int i=0; iGetNextCell ( numPoints, points ); + + vtkSmartPointer container = vtkSmartPointer::New(); + double* point = m_FiberPolyData->GetPoint(points[0]); + vtkIdType pointId = newPoints->InsertNextPoint(point); + container->GetPointIds()->InsertNextId(pointId); + float dtau = 0; + int cur_p = 1; + itk::Vector dR; + float normdR = 0; + for (;;) + { + while (dtau <= len && cur_p < numPoints) + { + itk::Vector v1; + point = m_FiberPolyData->GetPoint(points[cur_p-1]); + v1[0] = point[0]; + v1[1] = point[1]; + v1[2] = point[2]; + itk::Vector v2; + point = m_FiberPolyData->GetPoint(points[cur_p]); + v2[0] = point[0]; + v2[1] = point[1]; + v2[2] = point[2]; + + dR = v2 - v1; + normdR = std::sqrt(dR.GetSquaredNorm()); + dtau += normdR; + cur_p++; + } + + if (dtau >= len) + { + itk::Vector v1; + point = m_FiberPolyData->GetPoint(points[cur_p-1]); + v1[0] = point[0]; + v1[1] = point[1]; + v1[2] = point[2]; + + itk::Vector v2 = v1 - dR*( (dtau-len)/normdR ); + pointId = newPoints->InsertNextPoint(v2.GetDataPointer()); + container->GetPointIds()->InsertNextId(pointId); + } + else + { + point = m_FiberPolyData->GetPoint(points[numPoints-1]); + pointId = newPoints->InsertNextPoint(point); + container->GetPointIds()->InsertNextId(pointId); + break; + } + dtau = dtau-len; + } + + 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..0ee1a6e071 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(float len); + /*===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..72f1497f46 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->R = pVector(points[attrcnt*k]/spacing[0]-0.5, points[attrcnt*k+1]/spacing[1]-0.5,points[attrcnt*k+2]/spacing[2]-0.5); 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]; + index[1] = dp->R[1]; + index[2] = dp->R[2]; + 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/GibbsTracking/EnergyComputerBase.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/EnergyComputerBase.cpp index 9a2f4c5711..3b2d58cf82 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/EnergyComputerBase.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/EnergyComputerBase.cpp @@ -1,368 +1,369 @@ #ifndef _ENCOMPINTERFACE #define _ENCOMPINTERFACE #include "SphereInterpolator.cpp" #include "ParticleGrid.cpp" #include #include +#include inline float myATAN2(float y,float x) { float phi = acos(x); // float phi = ((x>=1.0) ? ((0.0000*x+-0.0000)) : ((x>=1.0) ? ((-10.0167*x+10.0167)) : ((x>=0.9) ? ((-3.1336*x+3.2713)) : ((x>=0.8) ? ((-1.9247*x+2.1833)) : ((x>=0.5) ? ((-1.3457*x+1.7200)) : ((x>=0.0) ? ((-1.0472*x+1.5708)) : ((x>=-0.5) ? ((-1.0472*x+1.5708)) : ((x>=-0.8) ? ((-1.3457*x+1.4216)) : ((x>=-0.9) ? ((-1.9247*x+0.9583)) : ((x>=-1.0) ? ((-3.1336*x+-0.1297)) : ((x>=-1.0) ? ((-10.0167*x+-6.8751)) : 1 ))))))))))); if (y < 0) phi = 2*PI - phi; if (phi<0) phi = phi + PI; return phi; } class EnergyComputerBase { public: float *m_QBallImageData; const int *m_QBallImageSize; SphereInterpolator *m_SphereInterpolator; ParticleGrid *m_ParticleGrid; int w,h,d; float voxsize_w; float voxsize_h; float voxsize_d; int w_sp,h_sp,d_sp; float voxsize_sp_w; float voxsize_sp_h; float voxsize_sp_d; int nip; // number of data vertices on sphere float *m_MaskImageData; float *cumulspatprob; int *spatidx; int scnt; float eigen_energy; vnl_matrix_fixed m_RotationMatrix; EnergyComputerBase(float *qBallImageData, const int *qBallImageSize, double *voxsize, SphereInterpolator *sp, ParticleGrid *pcon, float *maskImageData, int spmult, vnl_matrix_fixed rotMatrix) { m_RotationMatrix = rotMatrix; m_QBallImageData = qBallImageData; m_QBallImageSize = qBallImageSize; m_SphereInterpolator = sp; m_MaskImageData = maskImageData; nip = m_QBallImageSize[0]; w = m_QBallImageSize[1]; h = m_QBallImageSize[2]; d = m_QBallImageSize[3]; voxsize_w = voxsize[0]; voxsize_h = voxsize[1]; voxsize_d = voxsize[2]; w_sp = m_QBallImageSize[1]*spmult; h_sp = m_QBallImageSize[2]*spmult; d_sp = m_QBallImageSize[3]*spmult; voxsize_sp_w = voxsize[0]/spmult; voxsize_sp_h = voxsize[1]/spmult; voxsize_sp_d = voxsize[2]/spmult; fprintf(stderr,"Data size (voxels) : %i x %i x %i\n",w,h,d); fprintf(stderr,"voxel size: %f x %f x %f\n",voxsize_w,voxsize_h,voxsize_d); fprintf(stderr,"mask_oversamp_mult: %i\n",spmult); if (nip != sp->nverts) { fprintf(stderr,"EnergyComputer: error during init: data does not match with interpolation scheme\n"); } m_ParticleGrid = pcon; int totsz = w_sp*h_sp*d_sp; cumulspatprob = (float*) malloc(sizeof(float) * totsz); spatidx = (int*) malloc(sizeof(int) * totsz); if (cumulspatprob == 0 || spatidx == 0) { fprintf(stderr,"EnergyCOmputerBase: out of memory!\n"); return ; } scnt = 0; cumulspatprob[0] = 0; for (int x = 1; x < w_sp;x++) for (int y = 1; y < h_sp;y++) for (int z = 1; z < d_sp;z++) { int idx = x+(y+z*h_sp)*w_sp; if (m_MaskImageData[idx] > 0.5) { cumulspatprob[scnt+1] = cumulspatprob[scnt] + m_MaskImageData[idx]; spatidx[scnt] = idx; scnt++; } } for (int k = 0; k < scnt; k++) { cumulspatprob[k] /= cumulspatprob[scnt]; } fprintf(stderr,"#active voxels: %i (in mask units) \n",scnt); } ~EnergyComputerBase() { free(cumulspatprob); free(spatidx); } virtual void setParameters() { } void drawSpatPosition(pVector *R) { float r = mtrand.frand(); int j; int rl = 1; int rh = scnt; while(rh != rl) { j = rl + (rh-rl)/2; if (r < cumulspatprob[j]) { rh = j; continue; } if (r > cumulspatprob[j]) { rl = j+1; continue; } break; } R->SetXYZ(voxsize_sp_w*((float)(spatidx[rh-1] % w_sp) + mtrand.frand()), voxsize_sp_h*((float)((spatidx[rh-1]/w_sp) % h_sp) + mtrand.frand()), voxsize_sp_d*((float)(spatidx[rh-1]/(w_sp*h_sp)) + mtrand.frand())); } float SpatProb(pVector R) { int rx = int(R.GetX()/voxsize_sp_w); int ry = int(R.GetY()/voxsize_sp_h); int rz = int(R.GetZ()/voxsize_sp_d); if (rx >= 0 && rx < w_sp && ry >= 0 && ry < h_sp && rz >= 0 && rz < d_sp){ return m_MaskImageData[rx + w_sp* (ry + h_sp*rz)]; } else return 0; } /* inline float evaluateODF(pVector &R, pVector &N, float &len) { const int CU = 10; pVector Rs; float Dn = 0; int xint,yint,zint,spatindex; sinterp->getInterpolation(N); for (int i=-CU; i < CU;i++) { Rs = R + (N * len) * ((float)i/CU); xint = int(Rs.x); yint = int(Rs.y); zint = int(Rs.z); if (xint > 0 && xint < w-1 && yint > 0 && yint < h-1 && zint > 0 && zint < d-1) { spatindex = (xint + w*(yint+h*zint)) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2]); } } Dn /= (float)(2*CU+1); return Dn; } */ inline float evaluateODF(pVector &R, pVector &N, float &len) { const int CU = 10; pVector Rs; float Dn = 0; int xint,yint,zint,spatindex; vnl_vector_fixed n; n[0] = N[0]; n[1] = N[1]; n[2] = N[2]; n = m_RotationMatrix*n; m_SphereInterpolator->getInterpolation(n); for (int i=-CU; i <= CU;i++) { Rs = R + (N * len) * ((float)i/CU); float Rx = Rs[0]/voxsize_w-0.5; float Ry = Rs[1]/voxsize_h-0.5; float Rz = Rs[2]/voxsize_d-0.5; xint = int(floor(Rx)); yint = int(floor(Ry)); zint = int(floor(Rz)); if (xint >= 0 && xint < w-1 && yint >= 0 && yint < h-1 && zint >= 0 && zint < d-1) { float xfrac = Rx-xint; float yfrac = Ry-yint; float zfrac = Rz-zint; float weight; weight = (1-xfrac)*(1-yfrac)*(1-zfrac); spatindex = (xint + w*(yint+h*zint)) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(1-zfrac); spatindex = (xint+1 + w*(yint+h*zint)) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(1-zfrac); spatindex = (xint + w*(yint+1+h*zint)) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(1-yfrac)*(zfrac); spatindex = (xint + w*(yint+h*(zint+1))) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(yfrac)*(1-zfrac); spatindex = (xint+1 + w*(yint+1+h*zint)) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(zfrac); spatindex = (xint + w*(yint+1+h*(zint+1))) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(zfrac); spatindex = (xint+1 + w*(yint+h*(zint+1))) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(yfrac)*(zfrac); spatindex = (xint+1 + w*(yint+1+h*(zint+1))) *nip; Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; } } Dn *= 1.0/(2*CU+1); return Dn; } /* inline float evaluateODF(pVector &R, pVector &N, float &len) { R.storeXYZ(); float Rx = pVector::store[0]/voxsize_w; float Ry = pVector::store[1]/voxsize_h; float Rz = pVector::store[2]/voxsize_d; int xint = int(Rx); int yint = int(Ry); int zint = int(Rz); if (xint >= 0 && xint < w-1 && yint >= 0 && yint < h-1 && zint >= 0 && zint < d-1) { float xfrac = Rx-xint; float yfrac = Ry-yint; float zfrac = Rz-zint; sinterp->getInterpolation(N); float weight; int spatindex; float Dn = 0; weight = (1-xfrac)*(1-yfrac)*(1-zfrac); spatindex = (xint + w*(yint+h*zint)) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(1-zfrac); spatindex = (xint+1 + w*(yint+h*zint)) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(1-zfrac); spatindex = (xint + w*(yint+1+h*zint)) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; weight = (1-xfrac)*(1-yfrac)*(zfrac); spatindex = (xint + w*(yint+h*(zint+1))) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; weight = (xfrac)*(yfrac)*(1-zfrac); spatindex = (xint+1 + w*(yint+1+h*zint)) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(zfrac); spatindex = (xint + w*(yint+1+h*(zint+1))) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(zfrac); spatindex = (xint+1 + w*(yint+h*(zint+1))) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; weight = (xfrac)*(yfrac)*(zfrac); spatindex = (xint+1 + w*(yint+1+h*(zint+1))) *nip; Dn += (dataimg[spatindex + sinterp->idx[0]]*sinterp->interpw[0] + dataimg[spatindex + sinterp->idx[1]]*sinterp->interpw[1] + dataimg[spatindex + sinterp->idx[2]]* sinterp->interpw[2])*weight; return Dn; } return 0; } */ virtual inline float computeExternalEnergy(pVector &R, pVector &N, float &cap, float &len, Particle *dp) { return 0;} virtual inline float computeInternalEnergy(Particle *p1) {return 0;} virtual inline float computeInternalEnergyConnection(Particle *p1,int ep1) {return 0;} virtual inline float computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) {return 0;} }; #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMC_randshift.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMC_randshift.cpp index 0b83f23287..299e787783 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMC_randshift.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMC_randshift.cpp @@ -1,618 +1,620 @@ #include "ParticleGrid.cpp" #include "RJMCMCBase.cpp" #include class RJMCMC : public RJMCMCBase { public: float T_in ; float T_ex ; float dens; float p_birth; float p_death; float p_shift; float p_shiftopt; float p_cap; float p_con; float sigma_g; float gamma_g; float Z_g; float dthres; float nthres; float T_prop; float stopprobability; float del_prob; float len_def; float len_sig; float cap_def; float cap_sig; float externalEnergy; float internalEnergy; float m_ChempotParticle; Track TrackProposal, TrackBackup; SimpSamp simpsamp; RJMCMC(float *points,int numPoints, float *dimg, const int *dsz, double *voxsz, double cellsz) : RJMCMCBase(points,numPoints,dimg,dsz,voxsz,cellsz) { externalEnergy = 0; internalEnergy = 0; } void SetParameters(float Temp, int numit, float plen, float curv_hardthres, float chempot_particle) { m_Iterations = numit; p_birth = 0.25; p_death = 0.05; p_shift = 0.15; p_shiftopt = 0.1; p_con = 0.45; p_cap = 0.0; m_ChempotParticle = chempot_particle; float sum = p_birth+p_death+p_shift+p_shiftopt+p_con; p_birth /= sum; p_death /= sum; p_shift /= sum; p_shiftopt /= sum; T_in = Temp; T_ex = 0.01; dens = exp(-chempot_particle/T_in); len_def = plen; len_sig = 0.0; cap_def = 1.0; cap_sig = 0.0; // shift proposal sigma_g = len_def/8.0; gamma_g = 1/(sigma_g*sigma_g*2); Z_g = pow(2*PI*sigma_g,3.0/2.0)*(PI*sigma_g/len_def); // conn proposal dthres = len_def; nthres = curv_hardthres; T_prop = 0.5; dthres *= dthres; stopprobability = exp(-1/T_prop); del_prob = 0.1; } void SetTemperature(float temp) { T_in = temp; dens = exp(-m_ChempotParticle/T_in); } void IterateOneStep() { float randnum = mtrand.frand(); //randnum = 0; /////////////////////////////////////////////////////////////// //////// Birth Proposal /////////////////////////////////////////////////////////////// if (randnum < p_birth) { #ifdef TIMING tic(&birthproposal_time); birthstats.propose(); #endif pVector R; enc->drawSpatPosition(&R); + + //fprintf(stderr,"drawn: %f, %f, %f\n",R[0],R[1],R[2]); //R.setXYZ(20.5*3, 35.5*3, 1.5*3); pVector N; N.rand_sphere(); //N.setXYZ(1,0,0); float cap = cap_def - cap_sig*mtrand.frand(); float len = len_def;// + len_sig*(mtrand.frand()-0.5); Particle prop; prop.R = R; prop.N = N; prop.cap = cap; prop.len = len; float prob = dens * p_death /((p_birth)*(m_ParticleGrid.pcnt+1)); float ex_energy = enc->computeExternalEnergy(R,N,cap,len,0); float in_energy = enc->computeInternalEnergy(&prop); prob *= exp((in_energy/T_in+ex_energy/T_ex)) ; if (prob > 1 || mtrand.frand() < prob) { Particle *p = m_ParticleGrid.newParticle(R); if (p!=0) { p->R = R; p->N = N; p->cap = cap; p->len = len; #ifdef TIMING birthstats.accepted(); #endif m_AcceptedProposals++; } } #ifdef TIMING toc(&birthproposal_time); #endif } /////////////////////////////////////////////////////////////// //////// Death Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death) { if (m_ParticleGrid.pcnt > 0) { #ifdef TIMING tic(&deathproposal_time); deathstats.propose(); #endif int pnum = rand()%m_ParticleGrid.pcnt; Particle *dp = &(m_ParticleGrid.particles[pnum]); if (dp->pID == -1 && dp->mID == -1) { float ex_energy = enc->computeExternalEnergy(dp->R,dp->N,dp->cap,dp->len,dp); float in_energy = enc->computeInternalEnergy(dp); float prob = m_ParticleGrid.pcnt * (p_birth) /(dens*p_death); //*SpatProb(dp->R); prob *= exp(-(in_energy/T_in+ex_energy/T_ex)) ; if (prob > 1 || mtrand.frand() < prob) { m_ParticleGrid.remove(pnum); #ifdef TIMING deathstats.accepted(); #endif m_AcceptedProposals++; } } #ifdef TIMING toc(&deathproposal_time); #endif } } /////////////////////////////////////////////////////////////// //////// Cap change Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death+p_cap) { if (m_ParticleGrid.pcnt > 0) { int pnum = rand()%m_ParticleGrid.pcnt; Particle *p = &(m_ParticleGrid.particles[pnum]); Particle prop_p = *p; prop_p.cap = cap_def - cap_sig*mtrand.frand(); float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,prop_p.cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); //float in_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) // - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float prob = exp(ex_energy/T_ex); // * SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand.frand() < prob) { p->cap = prop_p.cap; m_AcceptedProposals++; } } } /////////////////////////////////////////////////////////////// //////// Shift Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death+p_shift+p_cap) { float energy = 0; if (m_ParticleGrid.pcnt > 0) { #ifdef TIMING tic(&shiftproposal_time); shiftstats.propose(); #endif int pnum = rand()%m_ParticleGrid.pcnt; Particle *p = &(m_ParticleGrid.particles[pnum]); Particle prop_p = *p; prop_p.R.distortn(sigma_g); prop_p.N.distortn(sigma_g/(2*p->len)); prop_p.N.normalize(); float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float in_energy = enc->computeInternalEnergy(&prop_p) - enc->computeInternalEnergy(p); float prob = exp(ex_energy/T_ex+in_energy/T_in); // * SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand.frand() < prob) { pVector Rtmp = p->R; pVector Ntmp = p->N; p->R = prop_p.R; p->N = prop_p.N; if (!m_ParticleGrid.tryUpdateGrid(pnum)) { p->R = Rtmp; p->N = Ntmp; } #ifdef TIMING shiftstats.accepted(); #endif m_AcceptedProposals++; } #ifdef TIMING toc(&shiftproposal_time); #endif } } else if (randnum < p_birth+p_death+p_shift+p_shiftopt+p_cap) { float energy = 0; if (m_ParticleGrid.pcnt > 0) { int pnum = rand()%m_ParticleGrid.pcnt; Particle *p = &(m_ParticleGrid.particles[pnum]); bool no_proposal = false; Particle prop_p = *p; if (p->pID != -1 && p->mID != -1) { Particle *plus = m_ParticleGrid.ID_2_address[p->pID]; int ep_plus = (plus->pID == p->ID)? 1 : -1; Particle *minus = m_ParticleGrid.ID_2_address[p->mID]; int ep_minus = (minus->pID == p->ID)? 1 : -1; prop_p.R = (plus->R + plus->N * (plus->len * ep_plus) + minus->R + minus->N * (minus->len * ep_minus))*0.5; prop_p.N = plus->R - minus->R; prop_p.N.normalize(); } else if (p->pID != -1) { Particle *plus = m_ParticleGrid.ID_2_address[p->pID]; int ep_plus = (plus->pID == p->ID)? 1 : -1; prop_p.R = plus->R + plus->N * (plus->len * ep_plus * 2); prop_p.N = plus->N; } else if (p->mID != -1) { Particle *minus = m_ParticleGrid.ID_2_address[p->mID]; int ep_minus = (minus->pID == p->ID)? 1 : -1; prop_p.R = minus->R + minus->N * (minus->len * ep_minus * 2); prop_p.N = minus->N; } else no_proposal = true; if (!no_proposal) { float cos = prop_p.N*p->N; float p_rev = exp(-((prop_p.R-p->R).norm_square() + (1-cos*cos))*gamma_g)/Z_g; float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float in_energy = enc->computeInternalEnergy(&prop_p) - enc->computeInternalEnergy(p); float prob = exp(ex_energy/T_ex+in_energy/T_in)*p_shift*p_rev/(p_shiftopt+p_shift*p_rev); //* SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand.frand() < prob) { pVector Rtmp = p->R; pVector Ntmp = p->N; p->R = prop_p.R; p->N = prop_p.N; if (!m_ParticleGrid.tryUpdateGrid(pnum)) { p->R = Rtmp; p->N = Ntmp; } m_AcceptedProposals++; } } } } else { if (m_ParticleGrid.pcnt > 0) { #ifdef TIMING tic(&connproposal_time); connstats.propose(); #endif int pnum = rand()%m_ParticleGrid.pcnt; Particle *p = &(m_ParticleGrid.particles[pnum]); EndPoint P; P.p = p; P.ep = (mtrand.frand() > 0.5)? 1 : -1; RemoveAndSaveTrack(P); if (TrackBackup.proposal_probability != 0) { MakeTrackProposal(P); float prob = (TrackProposal.energy-TrackBackup.energy)/T_in ; // prob = exp(prob)*(TrackBackup.proposal_probability) // /(TrackProposal.proposal_probability); prob = exp(prob)*(TrackBackup.proposal_probability * pow(del_prob,TrackProposal.length)) /(TrackProposal.proposal_probability * pow(del_prob,TrackBackup.length)); if (mtrand.frand() < prob) { ImplementTrack(TrackProposal); #ifdef TIMING connstats.accepted(); #endif m_AcceptedProposals++; } else { ImplementTrack(TrackBackup); } } else ImplementTrack(TrackBackup); #ifdef TIMING toc(&connproposal_time); #endif } } } void ImplementTrack(Track &T) { for (int k = 1; k < T.length;k++) { m_ParticleGrid.createConnection(T.track[k-1].p,T.track[k-1].ep,T.track[k].p,-T.track[k].ep); } } void RemoveAndSaveTrack(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; TrackBackup.track[cnt] = Current; EndPoint Next; for (;;) { Next.p = 0; if (Current.ep == 1) { if (Current.p->pID != -1) { Next.p = m_ParticleGrid.ID_2_address[Current.p->pID]; Current.p->pID = -1; m_ParticleGrid.concnt--; } } else if (Current.ep == -1) { if (Current.p->mID != -1) { Next.p = m_ParticleGrid.ID_2_address[Current.p->mID]; Current.p->mID = -1; m_ParticleGrid.concnt--; } } else { fprintf(stderr,"RJMCMC_randshift: Connection inconsistent 3\n"); break; } if (Next.p == 0) // no successor { Next.ep = 0; // mark as empty successor break; } else { if (Next.p->pID == Current.p->ID) { Next.p->pID = -1; Next.ep = 1; } else if (Next.p->mID == Current.p->ID) { Next.p->mID = -1; Next.ep = -1; } else { fprintf(stderr,"RJMCMC_randshift: Connection inconsistent 4\n"); break; } } ComputeEndPointProposalDistribution(Current); AccumProb *= (simpsamp.probFor(Next)); if (Next.p == 0) // no successor -> break break; energy += enc->computeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); Current = Next; Current.ep *= -1; cnt++; TrackBackup.track[cnt] = Current; if (mtrand.rand() > del_prob) { break; } } TrackBackup.energy = energy; TrackBackup.proposal_probability = AccumProb; TrackBackup.length = cnt+1; } void MakeTrackProposal(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; TrackProposal.track[cnt++] = Current; Current.p->label = 1; for (;;) { // next candidate is already connected if ((Current.ep == 1 && Current.p->pID != -1) || (Current.ep == -1 && Current.p->mID != -1)) break; // track too long if (cnt > 250) break; ComputeEndPointProposalDistribution(Current); // // no candidates anymore // if (simpsamp.isempty()) // break; int k = simpsamp.draw(); // stop tracking proposed if (k==0) break; EndPoint Next = simpsamp.objs[k]; float probability = simpsamp.probFor(k); // accumulate energy and proposal distribution energy += enc->computeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); AccumProb *= probability; // track to next endpoint Current = Next; Current.ep *= -1; Current.p->label = 1; // put label to avoid loops TrackProposal.track[cnt++] = Current; } TrackProposal.energy = energy; TrackProposal.proposal_probability = AccumProb; TrackProposal.length = cnt; // clear labels for (int j = 0; j < TrackProposal.length;j++) { TrackProposal.track[j].p->label = 0; } } void ComputeEndPointProposalDistribution(EndPoint P) { Particle *p = P.p; int ep = P.ep; float dist,dot; pVector R = p->R + (p->N * ep*p->len); m_ParticleGrid.computeNeighbors(R); simpsamp.clear(); simpsamp.add(stopprobability,EndPoint(0,0)); for (;;) { Particle *p2 = m_ParticleGrid.getNextNeighbor(); if (p2 == 0) break; if (p!=p2 && p2->label == 0) { if (p2->mID == -1) { dist = (p2->R - p2->N * p2->len - R).norm_square(); if (dist < dthres) { dot = p2->N*p->N * ep; if (dot > nthres) { float en = enc->computeInternalEnergyConnection(p,ep,p2,-1); simpsamp.add(exp(en/T_prop),EndPoint(p2,-1)); } } } if (p2->pID == -1) { dist = (p2->R + p2->N * p2->len - R).norm_square(); if (dist < dthres) { dot = p2->N*p->N * (-ep); if (dot > nthres) { float en = enc->computeInternalEnergyConnection(p,ep,p2,+1); simpsamp.add(exp(en/T_prop),EndPoint(p2,+1)); } } } } } } }; 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