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 67553574ed..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,775 +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_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; // 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; } } 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; typedef std::vector< itk::Point > FiberTractType; typedef std::vector< FiberTractType > FiberBundleType; vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); 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/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp index 2e679ea8f2..837669f448 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.cpp @@ -1,471 +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() { } /* * set computed fibers from tractography algorithms */ void mitk::FiberBundleX::SetFiberPolyData(vtkSmartPointer fiberPD) { if (fiberPD == NULL) this->m_FiberPolyData = vtkSmartPointer::New(); else this->m_FiberPolyData = fiberPD; m_isModified = true; } /* * return fiberbundle as vtkPolyData * Depending on processing of input fibers, this method returns * the latest processed fibers. */ vtkSmartPointer mitk::FiberBundleX::GetFiberPolyData() { return m_FiberPolyData; } /*=================================== *++++ PROCESSING WITH FIBERS +++++++ ====================================*/ void mitk::FiberBundleX::DoColorCodingOrientationbased() { //===== FOR WRITING A TEST ======================== // colorT size == tupelComponents * tupelElements // compare color results // to cover this code 100% also polydata needed, where colorarray already exists // + one fiber with exactly 1 point // + one fiber with 0 points //================================================= /* make sure that processing colorcoding is only called when necessary */ if ( m_FiberPolyData->GetPointData()->HasArray(COLORCODING_ORIENTATION_BASED) && m_FiberPolyData->GetNumberOfPoints() == m_FiberPolyData->GetPointData()->GetArray(COLORCODING_ORIENTATION_BASED)->GetNumberOfTuples() ) { // fiberstructure is already colorcoded MITK_INFO << " NO NEED TO REGENERATE COLORCODING! " ; return; } /* Finally, execute color calculation */ vtkPoints* extrPoints = m_FiberPolyData->GetPoints(); int numOfPoints = extrPoints->GetNumberOfPoints(); //colors and alpha value for each single point, RGBA = 4 components unsigned char rgba[4] = {0,0,0,0}; int componentSize = sizeof(rgba); vtkUnsignedCharArray * colorsT = vtkUnsignedCharArray::New(); colorsT->Allocate(numOfPoints * componentSize); colorsT->SetNumberOfComponents(componentSize); colorsT->SetName(COLORCODING_ORIENTATION_BASED); /* checkpoint: does polydata contain any fibers */ int numOfFibers = m_FiberPolyData->GetNumberOfLines(); if (numOfFibers < 1) { MITK_INFO << "\n ========= Number of Fibers is 0 and below ========= \n"; return; } /* extract single fibers of fiberBundle */ vtkCellArray* fiberList = m_FiberPolyData->GetLines(); fiberList->InitTraversal(); for (int fi=0; fiGetNextCell(pointsPerFiber, idList); // MITK_INFO << "Fib#: " << fi << " of " << numOfFibers << " pnts in fiber: " << pointsPerFiber ; /* single fiber checkpoints: is number of points valid */ if (pointsPerFiber > 1) { /* operate on points of single fiber */ for (int i=0; i 0) { /* The color value of the current point is influenced by the previous point and next point. */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; vnl_vector_fixed< double, 3 > diff; diff = (diff1 - diff2) / 2.0; diff.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff[2])); rgba[3] = (unsigned char) (255.0); } else if (i==0) { /* First point has no previous point, therefore only diff1 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > nextPntvtk(extrPoints->GetPoint(idList[i+1])[0], extrPoints->GetPoint(idList[i+1])[1], extrPoints->GetPoint(idList[i+1])[2]); vnl_vector_fixed< double, 3 > diff1; diff1 = currentPntvtk - nextPntvtk; diff1.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff1[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff1[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff1[2])); rgba[3] = (unsigned char) (255.0); } else if (i==pointsPerFiber-1) { /* Last point has no next point, therefore only diff2 is taken */ vnl_vector_fixed< double, 3 > currentPntvtk(extrPoints->GetPoint(idList[i])[0], extrPoints->GetPoint(idList[i])[1],extrPoints->GetPoint(idList[i])[2]); vnl_vector_fixed< double, 3 > prevPntvtk(extrPoints->GetPoint(idList[i-1])[0], extrPoints->GetPoint(idList[i-1])[1], extrPoints->GetPoint(idList[i-1])[2]); vnl_vector_fixed< double, 3 > diff2; diff2 = currentPntvtk - prevPntvtk; diff2.normalize(); rgba[0] = (unsigned char) (255.0 * std::abs(diff2[0])); rgba[1] = (unsigned char) (255.0 * std::abs(diff2[1])); rgba[2] = (unsigned char) (255.0 * std::abs(diff2[2])); rgba[3] = (unsigned char) (255.0); } colorsT->InsertTupleValue(idList[i], rgba); } //end for loop } else if (pointsPerFiber == 1) { /* a single point does not define a fiber (use vertex mechanisms instead */ continue; // colorsT->InsertTupleValue(0, rgba); } else { MITK_INFO << "Fiber with 0 points detected... please check your tractography algorithm!" ; continue; } }//end for loop m_FiberPolyData->GetPointData()->AddArray(colorsT); /*========================= - this is more relevant for renderer than for fiberbundleX datastructure - think about sourcing this to a explicit method which coordinates colorcoding */ this->SetColorCoding(COLORCODING_ORIENTATION_BASED); m_isModified = true; // =========================== //mini test, shall be ported to MITK TESTINGS! if (colorsT->GetSize() != numOfPoints*componentSize) { MITK_INFO << "ALLOCATION ERROR IN INITIATING COLOR ARRAY"; } } void mitk::FiberBundleX::DoGenerateFiberIds() { if (m_FiberPolyData == NULL) return; // for (int i=0; i<10000000; ++i) // { // if(i%500 == 0) // MITK_INFO << i; // } // MITK_INFO << "Generating Fiber Ids"; vtkSmartPointer idFiberFilter = vtkSmartPointer::New(); idFiberFilter->SetInput(m_FiberPolyData); idFiberFilter->CellIdsOn(); // idFiberFilter->PointIdsOn(); // point id's are not needed idFiberFilter->SetIdsArrayName(FIBER_ID_ARRAY); idFiberFilter->FieldDataOn(); idFiberFilter->Update(); m_FiberIdDataSet = idFiberFilter->GetOutput(); MITK_INFO << "Generating Fiber Ids...[done] | " << m_FiberIdDataSet->GetNumberOfCells(); } void mitk::FiberBundleX::UpdateFiberGeometry() { float min = itk::NumericTraits::min(); float max = itk::NumericTraits::max(); float b[] = {max, min, max, min, max, min}; vtkCellArray* cells = m_FiberPolyData->GetLines(); cells->InitTraversal(); for (int i=0; iGetNumberOfCells(); i++) { vtkCell* cell = m_FiberPolyData->GetCell(i); int p = cell->GetNumberOfPoints(); vtkPoints* points = cell->GetPoints(); for (int j=0; jGetPoint(j, p); if (p[0]b[1]) b[1]=p[0]; if (p[1]b[3]) b[3]=p[1]; if (p[2]b[5]) b[5]=p[2]; } } // provide some buffer space at borders for(int i=0; i<=4; i+=2){ b[i] -=10; } for(int i=1; i<=5; i+=2){ b[i] +=10; } mitk::Geometry3D::Pointer geometry = mitk::Geometry3D::New(); geometry->SetImageGeometry(true); geometry->SetFloatBounds(b); this->SetGeometry(geometry); } /*============================== *++++ FIBER INFORMATION +++++++ ===============================*/ QStringList mitk::FiberBundleX::GetAvailableColorCodings() { QStringList availableColorCodings; int numColors = m_FiberPolyData->GetPointData()->GetNumberOfArrays(); for(int i=0; iGetPointData()->GetArrayName(i)); } //this controlstructure shall be implemented by the calling method if (availableColorCodings.isEmpty()) MITK_INFO << "no colorcodings available in fiberbundleX"; // for(int i=0; im_currentColorCoding; } void mitk::FiberBundleX::SetColorCoding(const char* requestedColorCoding) { // MITK_INFO << "FbX try to set colorCoding: " << requestedColorCoding << " compare with: " << COLORCODING_ORIENTATION_BASED; if(strcmp (COLORCODING_ORIENTATION_BASED,requestedColorCoding) == 0 ) { this->m_currentColorCoding = (char*) COLORCODING_ORIENTATION_BASED; this->m_isModified = true; } else if(strcmp (COLORCODING_FA_BASED,requestedColorCoding) == 0 ) { this->m_currentColorCoding = (char*) COLORCODING_FA_BASED; this->m_isModified = true; } else { MITK_INFO << "FIBERBUNDLE X: UNKNOWN COLORCODING in FIBERBUNDLEX Datastructure"; this->m_currentColorCoding = "---"; //will cause blank colorcoding of fibers this->m_isModified = true; } } bool mitk::FiberBundleX::isFiberBundleXModified() { return m_isModified; } void mitk::FiberBundleX::setFBXModificationDone() { m_isModified = false; } -void mitk::FiberBundleX::ResampleFibers() +// 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 cellArray = m_FiberPolyData->GetLines(); - for (int i=0; iGetNumberOfLines(); i++) + vtkSmartPointer vLines = m_FiberPolyData->GetLines(); + vLines->InitTraversal(); + int numberOfLines = vLines->GetNumberOfCells(); + + for (int i=0; iGetCell(cellId, numPoints, points); + vtkIdType numPoints(0); + vtkIdType* points(NULL); + vLines->GetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); - // insert start double* point = m_FiberPolyData->GetPoint(points[0]); vtkIdType pointId = newPoints->InsertNextPoint(point); container->GetPointIds()->InsertNextId(pointId); - float Leng = 0; float dtau = 0; int cur_p = 1; - int cur_i = 1; - pVector dR; - float normdR; + itk::Vector dR; + float normdR = 0; + for (;;) { while (dtau <= len && cur_p < numPoints) { - dR = points[cur_p] - points[cur_p-1]; - normdR = sqrt(dR.norm_square()); + 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; - Leng += normdR; cur_p++; } - if (dtau >= len) // if particles reach next voxel + if (dtau >= len) { - point = points[cur_p-1] - dR*( (dtau-len)/normdR ); - pointId = newPoints->InsertNextPoint(point); + 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; - cur_i++; } newCellArray->InsertNextCell(container); } newPoly->SetPoints(newPoints); newPoly->SetLines(newCellArray); m_FiberPolyData = newPoly; UpdateFiberGeometry(); } /* ESSENTIAL IMPLEMENTATION OF SUPERCLASS METHODS */ void mitk::FiberBundleX::UpdateOutputInformation() { } void mitk::FiberBundleX::SetRequestedRegionToLargestPossibleRegion() { } bool mitk::FiberBundleX::RequestedRegionIsOutsideOfTheBufferedRegion() { return false; } bool mitk::FiberBundleX::VerifyRequestedRegion() { return true; } void mitk::FiberBundleX::SetRequestedRegion( itk::DataObject *data ) { } diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h index 207d4c3fd3..0ee1a6e071 100644 --- a/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundleX/mitkFiberBundleX.h @@ -1,129 +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.) 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(); + 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 50b8b03c5f..72f1497f46 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp @@ -1,155 +1,155 @@ #ifndef _BUILDFIBRES #define _BUILDFIBRES #include #include #include #include using namespace std; #define PI 3.1415926536 #include "ParticleGrid.cpp" #include #include #include #include #include #include #include class FiberBuilder { public: Particle *particles; int pcnt; int attrcnt; typedef vector< Particle* > ParticleContainerType; typedef vector< ParticleContainerType* > FiberContainerType; vtkSmartPointer m_VtkCellArray; vtkSmartPointer m_VtkPoints; typedef itk::Vector OdfVectorType; typedef itk::Image ItkQBallImgType; ItkQBallImgType::Pointer m_ItkQBallImage; FiberBuilder(float *points, int numPoints, double spacing[], ItkQBallImgType::Pointer image) { 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(); } ~FiberBuilder() { free(particles); } 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) { vtkSmartPointer container = vtkSmartPointer::New(); dp->label = cur_label; dp->numerator = 0; labelPredecessors(dp, container); labelSuccessors(dp, container); cur_label++; if(container->GetPointIds()->GetNumberOfIds()>minSize) { m_VtkCellArray->InsertNextCell(container); numFibers++; } } } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(m_VtkPoints); fiberPolyData->SetLines(m_VtkCellArray); return fiberPolyData; } void AddPoint(Particle *dp, vtkSmartPointer container) { itk::ContinuousIndex index; - index[0] = dp->R[0]-0.5; - index[1] = dp->R[1]-0.5; - index[2] = dp->R[2]-0.5; + 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, 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); } } AddPoint(dp, container); } void labelSuccessors(Particle *dp, vtkSmartPointer container) { 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); } } } }; #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)); } } } } } } };