diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake index 1221a5b610..11ea4a277a 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake @@ -1,72 +1,76 @@ SET(SRC_CPP_FILES QmitkODFDetailsWidget.cpp QmitkODFRenderWidget.cpp ) SET(INTERNAL_CPP_FILES mitkPluginActivator.cpp QmitkQBallReconstructionView.cpp QmitkPreprocessingView.cpp QmitkDiffusionDicomImportView.cpp QmitkDiffusionQuantificationView.cpp QmitkTensorReconstructionView.cpp QmitkDiffusionImagingPublicPerspective.cpp QmitkControlVisualizationPropertiesView.cpp QmitkODFDetailsView.cpp + QmitkGlobalFiberTrackingView.cpp ) SET(UI_FILES src/internal/QmitkQBallReconstructionViewControls.ui src/internal/QmitkPreprocessingViewControls.ui src/internal/QmitkDiffusionDicomImportViewControls.ui src/internal/QmitkDiffusionQuantificationViewControls.ui src/internal/QmitkTensorReconstructionViewControls.ui src/internal/QmitkControlVisualizationPropertiesViewControls.ui src/internal/QmitkODFDetailsViewControls.ui + src/internal/QmitkGlobalFiberTrackingViewControls.ui ) SET(MOC_H_FILES src/internal/mitkPluginActivator.h src/internal/QmitkQBallReconstructionView.h src/internal/QmitkPreprocessingView.h src/internal/QmitkDiffusionDicomImportView.h src/internal/QmitkDiffusionImagingPublicPerspective.h src/internal/QmitkDiffusionQuantificationView.h src/internal/QmitkTensorReconstructionView.h src/internal/QmitkControlVisualizationPropertiesView.h src/internal/QmitkODFDetailsView.h src/QmitkODFRenderWidget.h src/QmitkODFDetailsWidget.h + src/internal/QmitkGlobalFiberTrackingView.h ) SET(CACHED_RESOURCE_FILES # list of resource files which can be used by the plug-in # system without loading the plug-ins shared library, # for example the icon used in the menu and tabs for the # plug-in views in the workbench plugin.xml resources/preprocessing.png resources/dwiimport.png resources/quantification.png resources/reconodf.png resources/recontensor.png resources/vizControls.png - resources/odf.png + resources/OdfDetails.png + resources/GlobalTracking.png ) SET(QRC_FILES # uncomment the following line if you want to use Qt resources resources/QmitkDiffusionImaging.qrc ) SET(CPP_FILES ) foreach(file ${SRC_CPP_FILES}) SET(CPP_FILES ${CPP_FILES} src/${file}) endforeach(file ${SRC_CPP_FILES}) foreach(file ${INTERNAL_CPP_FILES}) SET(CPP_FILES ${CPP_FILES} src/internal/${file}) endforeach(file ${INTERNAL_CPP_FILES}) diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml index e75f2c811b..fe30ab1d8a 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml @@ -1,61 +1,68 @@ + icon="resources/OdfDetails.png" /> + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/GlobalTracking.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/GlobalTracking.png new file mode 100644 index 0000000000..d875ea96bd Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/GlobalTracking.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/OdfDetails.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/OdfDetails.png new file mode 100644 index 0000000000..c1eff6a986 Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/OdfDetails.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingView.cpp new file mode 100644 index 0000000000..30e297870c --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingView.cpp @@ -0,0 +1,699 @@ +/*========================================================================= +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 "QmitkGlobalFiberTrackingView.h" +#include + +// Qt +#include +#include +#include + +// MITK +#include +#include +#include +#include + +// ITK +#include +#include + +// MISC +#include + + +QmitkTrackingWorker::QmitkTrackingWorker(QmitkGlobalFiberTrackingView* view) + : m_View(view) +{ + +} + +void QmitkTrackingWorker::run() +{ + 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(); + + m_View->m_GlobalTracker = QmitkGlobalFiberTrackingView::GlobalTrackingFilterType::New(); + m_View->m_GlobalTracker->SetInput0(m_View->m_ItkQBallImage.GetPointer()); + m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); + 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 QmitkGlobalFiberTrackingView::VIEW_ID = +"org.mitk.views.globalfibertracking"; + +QmitkGlobalFiberTrackingView::QmitkGlobalFiberTrackingView() + : QmitkFunctionality() + , m_Controls( 0 ) + , m_MultiWidget( NULL ) + , m_ThreadIsRunning(false) + , m_GlobalTracker(NULL) + , m_QBallImage(NULL) + , m_MaskImage(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); +} + +QmitkGlobalFiberTrackingView::~QmitkGlobalFiberTrackingView() +{ + delete m_TrackingTimer; +} + +// update tracking status and generate fiber bundle +void QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::StopGlobalTracking() +{ + if (m_GlobalTracker.IsNull()) + return; + m_GlobalTracker->SetAbortTracking(true); + m_Controls->m_GlobalTrackingStop->setEnabled(false); + m_Controls->m_GlobalTrackingStop->setText("Stopping Tractography ..."); +} + +// update gui elements and generate fiber bundle after tracking is finished +void QmitkGlobalFiberTrackingView::AfterThread() +{ + m_ThreadIsRunning = false; + m_TrackingTimer->stop(); + + UpdateGUI(); + UpdateTrackingStatus(); + GenerateFiberBundle(); +} + +// start tracking timer and update gui elements before tracking is started +void QmitkGlobalFiberTrackingView::BeforeThread() +{ + m_ThreadIsRunning = true; + m_TrackingTime = QTime::currentTime(); + m_ElapsedTime = 0; + m_TrackingTimer->start(1000); + + UpdateGUI(); +} + +// setup gui elements and signal/slot connections +void QmitkGlobalFiberTrackingView::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::QmitkGlobalFiberTrackingViewControls; + m_Controls->setupUi( parent ); + + AdvancedSettings(); + + connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); + connect( m_Controls->m_GlobalTrackingStop, SIGNAL(clicked()), this, SLOT(StopGlobalTracking()) ); + connect( m_Controls->m_GlobalTrackingStart, SIGNAL(clicked()), this, SLOT(StartGlobalTracking()) ); + connect( m_Controls->m_SetMaskButton, SIGNAL(clicked()), this, SLOT(SetMask()) ); + 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 QmitkGlobalFiberTrackingView::SetInExBalance(int value) +{ + m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); +} + +void QmitkGlobalFiberTrackingView::SetFiberLength(int value) +{ + m_Controls->m_FiberLengthLabel->setText(QString::number(value)); +} + +void QmitkGlobalFiberTrackingView::SetParticleWeight(int value) +{ + m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); +} + +void QmitkGlobalFiberTrackingView::SetStartTemp(int value) +{ + m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); +} + +void QmitkGlobalFiberTrackingView::SetEndTemp(int value) +{ + m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); +} + +void QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::StdMultiWidgetAvailable(QmitkStdMultiWidget &stdMultiWidget) +{ + m_MultiWidget = &stdMultiWidget; +} + +void QmitkGlobalFiberTrackingView::StdMultiWidgetNotAvailable() +{ + m_MultiWidget = NULL; +} + +// called if datamanager selection changes +void QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::UpdateGUI() +{ + if (!m_ThreadIsRunning && m_QBallSelected) + { + m_Controls->m_GlobalTrackingStop->setEnabled(false); + m_Controls->m_GlobalTrackingStart->setEnabled(true); + m_Controls->m_LoadTrackingParameters->setEnabled(true); + m_Controls->m_MaskFrame->setEnabled(true); + m_Controls->m_IterationsSlider->setEnabled(true); + m_Controls->m_AdvancedFrame->setEnabled(true); + m_Controls->m_GlobalTrackingStop->setText("Stop Tractography"); + m_Controls->m_GlobalTrackingStart->setToolTip("Start tractography. No further change of parameters possible."); + m_Controls->m_GlobalTrackingStop->setToolTip(""); + } + else if (!m_ThreadIsRunning) + { + m_Controls->m_GlobalTrackingStop->setEnabled(false); + m_Controls->m_GlobalTrackingStart->setEnabled(false); + m_Controls->m_LoadTrackingParameters->setEnabled(true); + m_Controls->m_MaskFrame->setEnabled(true); + m_Controls->m_IterationsSlider->setEnabled(true); + m_Controls->m_AdvancedFrame->setEnabled(true); + m_Controls->m_GlobalTrackingStop->setText("Stop Tractography"); + m_Controls->m_GlobalTrackingStart->setToolTip("No Q-Ball image selected."); + m_Controls->m_GlobalTrackingStop->setToolTip(""); + } + else + { + m_Controls->m_GlobalTrackingStop->setEnabled(true); + m_Controls->m_GlobalTrackingStart->setEnabled(false); + m_Controls->m_LoadTrackingParameters->setEnabled(false); + m_Controls->m_MaskFrame->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_GlobalTrackingStart->setToolTip("Tracking in progress."); + m_Controls->m_GlobalTrackingStop->setToolTip("Stop tracking and display results."); + } +} + +// show/hide advanced settings frame +void QmitkGlobalFiberTrackingView::AdvancedSettings() +{ + m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); +} + +// set mask image data node +void QmitkGlobalFiberTrackingView::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; + } + } +} + +// cast image to float +template +void QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::StartGlobalTracking() +{ + if(m_ThreadIsRunning) + { + MITK_WARN("QmitkGlobalFiberTrackingView")<<"Thread already running!"; + return; + } + + if (!m_QBallSelected) + { + // Nothing selected. Inform the user and return + QMessageBox::information( NULL, "Template", "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? + if(m_Controls->m_MaskImageEdit->text().compare("N/A") != 0) + { + m_MaskImage = 0; + mitk::BaseData* data = m_MaskImageNode->GetData(); + if (data) + { + // test if this data item is an image or not (could also be a surface or something totally different) + mitk::Image* tmpImage = dynamic_cast( data ); + if (tmpImage) + { + mitk::Image::Pointer mitkMaskImg = mitk::Image::New(); + AccessFixedDimensionByItk_1(tmpImage, CastToFloat, 3, mitkMaskImg); + typedef mitk::ImageToItk CastType; + CastType::Pointer caster = CastType::New(); + caster->SetInput(mitkMaskImg); + caster->Update(); + m_MaskImage = caster->GetOutput(); + } + } + } + + // 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); + } + } + + 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 QmitkGlobalFiberTrackingView::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(); + + 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()); + + 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 QmitkGlobalFiberTrackingView::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", m_Iterations); + 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 QmitkGlobalFiberTrackingView::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 QmitkGlobalFiberTrackingView::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/QmitkGlobalFiberTrackingView.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingView.h new file mode 100644 index 0000000000..6669acf79a --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingView.h @@ -0,0 +1,162 @@ +/*========================================================================= +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 QmitkGlobalFiberTrackingView_h +#define QmitkGlobalFiberTrackingView_h + +#include + +#include + +#include "ui_QmitkGlobalFiberTrackingViewControls.h" + +#include +#include +#include +#include +#include + +class QmitkGlobalFiberTrackingView; + +class QmitkTrackingWorker : public QObject +{ + Q_OBJECT + +public: + + QmitkTrackingWorker(QmitkGlobalFiberTrackingView* view); + +public slots: + + void run(); + +private: + + QmitkGlobalFiberTrackingView* m_View; +}; + +/*! + \brief QmitkGlobalFiberTrackingView + + \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 GlobalTractographyFilter; +} + +class QmitkGlobalFiberTrackingView : 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::GlobalTractographyFilter GlobalTrackingFilterType; + + static const std::string VIEW_ID; + + QmitkGlobalFiberTrackingView(); + virtual ~QmitkGlobalFiberTrackingView(); + + virtual void CreateQtPartControl(QWidget *parent); + + virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); + virtual void StdMultiWidgetNotAvailable(); + +signals: + +protected slots: + + void StartGlobalTracking(); + void StopGlobalTracking(); + void AfterThread(); + void BeforeThread(); + void TimerUpdate(); + void SetMask(); + 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::QmitkGlobalFiberTrackingViewControls* m_Controls; + QmitkStdMultiWidget* m_MultiWidget; + + // data objects + mitk::FiberBundle::Pointer m_FiberBundle; + MaskImgType::Pointer m_MaskImage; + 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_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 // _QMITKGLOBALFIBERTRACKINGVIEW_H_INCLUDED + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingViewControls.ui new file mode 100644 index 0000000000..7f907bfd2b --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGlobalFiberTrackingViewControls.ui @@ -0,0 +1,946 @@ + + + QmitkGlobalFiberTrackingViewControls + + + + 0 + 0 + 463 + 1011 + + + + + 0 + 0 + + + + + 0 + 0 + + + + QmitkTemplate + + + + 0 + + + 3 + + + 0 + + + 3 + + + 0 + + + + + + + + false + + + + 0 + + + 0 + + + + + true + + + false + + + QFrame::NoFrame + + + QFrame::Raised + + + + QFormLayout::AllNonFixedFieldsGrow + + + 0 + + + 4 + + + 0 + + + + + Mask Image + + + + + + + Mask image. Particles will be sampled with a probability according to the voxel value. + + + QFrame::StyledPanel + + + QFrame::Raised + + + + 0 + + + + + true + + + => + + + + + + + true + + + N/A + + + true + + + + + + + + + + + + + + + + + + + Iterations: 10^7 + + + + + + + true + + + Visualize Tractography + + + true + + + + + + + Advanced Settings + + + + + + + true + + + QFrame::StyledPanel + + + QFrame::Raised + + + + 0 + + + 4 + + + 0 + + + + + + + + + + + + + + auto + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + + + + + + auto = 1.5 * min. spacing; l + + + 100 + + + 1 + + + Qt::Horizontal + + + QSlider::NoTicks + + + + + + + + + + + + + + + + auto + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + + + + + + auto = 0.5 * min. spacing; sigma + + + 100 + + + 1 + + + Qt::Horizontal + + + QSlider::NoTicks + + + + + + + + + + + + + + + + 0.001 + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + + + + + + unitless w + + + 1 + + + 1000 + + + 1 + + + 10 + + + Qt::Horizontal + + + true + + + QSlider::NoTicks + + + + + + + 0.1 + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + + + + + + 1 + + + 100 + + + 1 + + + 10 + + + Qt::Horizontal + + + false + + + false + + + QSlider::NoTicks + + + + + + + 0.001 + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + + + + + + 1 + + + 99 + + + 1 + + + 10 + + + Qt::Horizontal + + + QSlider::NoTicks + + + + + + + + + + + + + + + + 0 + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + + + + + + IE Bias < 0 < EE Bias + + + -50 + + + 50 + + + 1 + + + Qt::Horizontal + + + QSlider::NoTicks + + + + + + + + + + + + + + + + 10 + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + + + + + + + Only fibers containing more than the specified number of particles are accepted. + + + 100 + + + 1 + + + 10 + + + Qt::Horizontal + + + QSlider::NoTicks + + + + + + + + + + + + + + + + Particle Length: + + + + + + + + + + + + + + + + Particle Width: + + + + + + + + + + + + + + + + Particle Weight: + + + + + + + Start Temperature: + + + + + + + End Temperature: + + + + + + + + + + + + + + + + Balance In/Ex Energy: + + + + + + + + + + + + + + + + Min. Fiber Length: + + + + + + + true + + + Use mean subtracted ODFs (recommended). + + + Subtract ODF Mean + + + true + + + + + + + Qt::Horizontal + + + QSizePolicy::Fixed + + + + 60 + 20 + + + + + + + + + + + true + + + Save current parameters as xml (.gtp) + + + Qt::LeftToRight + + + Save Parameters + + + + :/qmitk/btnMoveDown.png:/qmitk/btnMoveDown.png + + + + + + + true + + + Load parameters from xml file (.gtp) + + + Qt::LeftToRight + + + Load Parameters + + + + :/qmitk/btnMoveUp.png:/qmitk/btnMoveUp.png + + + + + + + Specify number of iterations for the tracking algorithm. + + + 11 + + + 6 + + + Qt::Horizontal + + + QSlider::TicksBelow + + + + + + + false + + + No Q-Ball image selected. + + + Qt::LeftToRight + + + Start Tractography + + + + :/qmitk/play.xpm:/qmitk/play.xpm + + + + + + + false + + + + + + Qt::LeftToRight + + + Stop Tractography + + + + :/qmitk/stop.xpm:/qmitk/stop.xpm + + + + + + + + + + + 0 + 0 + + + + + 0 + 0 + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + 0 + + + 7 + + + 0 + + + + + Will only be updated if tracking is visualized + + + Will only be updated if tracking is visualized + + + + + + Accepted Fibers: + + + + + + + Will only be updated if tracking is visualized + + + + + + + + + - + + + + + + + + + + + + + + + + Progress: + + + + + + + + + + + + + + + + - + + + + + + + + + + + + + + + + Connections: + + + + + + + + + + + + + + + + Particles: + + + + + + + + + + + + + + + + - + + + + + + + + + + + + + + + + - + + + + + + + + + + + + + + + + Tracking Time: + + + + + + + + + + + + + + + + - + + + + + + + + + + + + + + + + Proposal Acceptance Rate: + + + + + + + + + + + + + + + + - + + + + + + + + + + + + + Qt::Vertical + + + QSizePolicy::Expanding + + + + 0 + 0 + + + + + + + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp index 3413428379..1715b75d19 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp @@ -1,38 +1,40 @@ #include "mitkPluginActivator.h" #include #include "src/internal/QmitkDiffusionImagingPublicPerspective.h" #include "src/internal/QmitkQBallReconstructionView.h" #include "src/internal/QmitkPreprocessingView.h" #include "src/internal/QmitkDiffusionDicomImportView.h" #include "src/internal/QmitkDiffusionQuantificationView.h" #include "src/internal/QmitkTensorReconstructionView.h" #include "src/internal/QmitkControlVisualizationPropertiesView.h" #include "src/internal/QmitkODFDetailsView.h" +#include "src/internal/QmitkGlobalFiberTrackingView.h" namespace mitk { void PluginActivator::start(ctkPluginContext* context) { BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionImagingPublicPerspective, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkQBallReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkPreprocessingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionDicomImport, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionQuantificationView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkTensorReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkControlVisualizationPropertiesView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkODFDetailsView, context) + BERRY_REGISTER_EXTENSION_CLASS(QmitkGlobalFiberTrackingView, context) } void PluginActivator::stop(ctkPluginContext* context) { Q_UNUSED(context) } } Q_EXPORT_PLUGIN2(org_mitk_gui_qt_diffusionimaging, mitk::PluginActivator) diff --git a/Modules/DiffusionImaging/CMakeLists.txt b/Modules/DiffusionImaging/CMakeLists.txt index 63f18398f3..7766735ede 100644 --- a/Modules/DiffusionImaging/CMakeLists.txt +++ b/Modules/DiffusionImaging/CMakeLists.txt @@ -1,20 +1,26 @@ FIND_PACKAGE(ITK) IF(ITK_GDCM_DIR) INCLUDE(${ITK_GDCM_DIR}/GDCMConfig.cmake) IF(GDCM_MAJOR_VERSION EQUAL 2) ADD_DEFINITIONS(-DGDCM2) SET(ITK_USES_GDCM2 1) ENDIF(GDCM_MAJOR_VERSION EQUAL 2) ENDIF(ITK_GDCM_DIR) MITK_CREATE_MODULE( MitkDiffusionImaging SUBPROJECTS MITK-DTI INCLUDE_DIRS Algorithms DicomImport Interactions IODataStructures/DiffusionWeightedImages IODataStructures/QBallImages IODataStructures/TensorImages IODataStructures/FiberBundle IODataStructures/PlanarFigureComposite IODataStructures Reconstruction Tractography Rendering ${CMAKE_CURRENT_BINARY_DIR} DEPENDS MitkExt SceneSerializationBase QmitkExt PACKAGE_DEPENDS Boost ) +MITK_USE_MODULE(MitkDiffusionImaging) +if(MitkDiffusionImaging_IS_ENABLED) + file(DOWNLOAD http://mitk.org/download/data/FibertrackingLUT.tar.gz ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/FibertrackingLUT.tar.gz TIMEOUT 10) + execute_process(COMMAND cmake -E chdir ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} tar xzf FibertrackingLUT.tar.gz) +endif() + ADD_SUBDIRECTORY(Testing) - + CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/mitkDiffusionImagingConfigure.h.in ${CMAKE_CURRENT_BINARY_DIR}/mitkDiffusionImagingConfigure.h) diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticle.cpp b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticle.cpp new file mode 100644 index 0000000000..d2da8870ee --- /dev/null +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticle.cpp @@ -0,0 +1,45 @@ +/*========================================================================= + 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 "mitkParticle.h" + +namespace mitk +{ + Particle::Particle() + { + + } + + Particle::~Particle() + { + + } + + void Particle::SetDirection(VectorType dir) + { + this->m_Direction = dir; + } + + void Particle::SetPosition(VectorType pos) + { + this->m_Position = pos; + } + + Particle::VectorType Particle::GetDirection() + { + return this->m_Direction; + } + + Particle::VectorType Particle::GetPosition() + { + return this->m_Position; + } +} diff --git a/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticle.h b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticle.h new file mode 100644 index 0000000000..d8f5237138 --- /dev/null +++ b/Modules/DiffusionImaging/IODataStructures/FiberBundle/mitkParticle.h @@ -0,0 +1,51 @@ + +/*========================================================================= + 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. + + =========================================================================*/ + + +#ifndef _MITK_Particle_H +#define _MITK_Particle_H + +#include "mitkBaseData.h" +#include "MitkDiffusionImagingExports.h" + +namespace mitk { + + /** + * \brief Base Class for Fiber Segments; */ + class MitkDiffusionImaging_EXPORT Particle + { + public: + Particle(); + virtual ~Particle(); + + typedef vnl_vector_fixed VectorType; + + void SetDirection(VectorType dir); + void SetPosition(VectorType pos); + + VectorType GetPosition(); + VectorType GetDirection(); + + protected: + VectorType m_Position; // position + VectorType m_Direction; // direction + }; + +} // namespace mitk + +#endif /* _MITK_Particle_H */ diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/AccumulateBilin.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/AccumulateBilin.cpp new file mode 100644 index 0000000000..1f78d15d7b --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/AccumulateBilin.cpp @@ -0,0 +1,87 @@ +#include +#include "mex.h" +#include "matrix.h" +#define REAL float + + + + + + + + +void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) +{ + + if(nrhs!=2) { + mexPrintf("wrong usage!!\n",nrhs); + return; + } else if(nlhs>4) { + printf("Too many output arguments\n"); + return; + } + + int pcnt = 0; + const mxArray *Dim; + Dim = prhs[pcnt++]; + REAL *dim = (REAL*) mxGetData(Dim); + int w = (int) dim[0]; + int h = (int) dim[1]; + int d = (int) dim[2]; + + const mxArray *Pts; + Pts = prhs[pcnt++]; + const int numdim = mxGetNumberOfDimensions(Pts); + const int *pdims = mxGetDimensions(Pts); + int numPts = pdims[1]; + REAL *pts = (REAL*) mxGetData(Pts); + + + int dims[3]; + dims[0] = w; + dims[1] = h; + dims[2] = d; + plhs[0] = mxCreateNumericArray(3,dims,mxGetClassID(Dim),mxREAL); + REAL *accu = (REAL*) mxGetData(plhs[0]); + + + for (int i = 0; i < numPts; i++) + { + int idx = 3*i; + + int px = (int) (pts[idx]); + if (px < 0 || px >= w-1) + continue; + int py = (int) (pts[idx+1]); + if (py < 0 || py >= h-1) + continue; + int pz = (int) (pts[idx+2]); + if (pz < 0 || pz >= d-1) + continue; + float frac_x = pts[idx ] - px; + float frac_y = pts[idx+1] - py; + float frac_z = pts[idx+2] - pz; + + + accu[px + w*(py+h*pz)] += (1-frac_x)*(1-frac_y)*(1-frac_z); + + accu[px+1 + w*(py+h*pz)] += (frac_x)*(1-frac_y)*(1-frac_z); + accu[px + w*(py+1+h*pz)] += (1-frac_x)*(frac_y)*(1-frac_z); + accu[px + w*(py+h*pz+h)] += (1-frac_x)*(1-frac_y)*(frac_z); + + accu[px + w*(py+1+h*pz+h)] += (1-frac_x)*(frac_y)*(frac_z); + accu[px+1 + w*(py+h*pz+h)] += (frac_x)*(1-frac_y)*(frac_z); + accu[px+1 + w*(py+1+h*pz)] += (frac_x)*(frac_y)*(1-frac_z); + + accu[px+1 + w*(py+1+h*pz+h)] += (frac_x)*(frac_y)*(frac_z); + + } + + +} + + + + + + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/AccumulateBilinWeighted.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/AccumulateBilinWeighted.cpp new file mode 100644 index 0000000000..c60c7c9f6c --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/AccumulateBilinWeighted.cpp @@ -0,0 +1,94 @@ +#include +#include "mex.h" +#include "matrix.h" +#define REAL float + + + + + + + + +void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) +{ + + if(nrhs!=3) { + mexPrintf("wrong usage!!\n",nrhs); + return; + } else if(nlhs>4) { + printf("Too many output arguments\n"); + return; + } + + int pcnt = 0; + const mxArray *Dim; + Dim = prhs[pcnt++]; + REAL *dim = (REAL*) mxGetData(Dim); + int w = (int) dim[0]; + int h = (int) dim[1]; + int d = (int) dim[2]; + + const mxArray *Pts; + Pts = prhs[pcnt++]; + const int numdim = mxGetNumberOfDimensions(Pts); + const int *pdims = mxGetDimensions(Pts); + int numPts = pdims[1]; + REAL *pts = (REAL*) mxGetData(Pts); + + const mxArray *Weights; + Weights = prhs[pcnt++]; + const int wnumdim = mxGetNumberOfDimensions(Weights); + const int *wpdims = mxGetDimensions(Weights); + int numW = wpdims[0]; + REAL *weights = (REAL*) mxGetData(Weights); + + + int dims[3]; + dims[0] = w; + dims[1] = h; + dims[2] = d; + plhs[0] = mxCreateNumericArray(3,dims,mxGetClassID(Dim),mxREAL); + REAL *accu = (REAL*) mxGetData(plhs[0]); + + + for (int i = 0; i < numPts; i++) + { + int idx = 3*i; + + int px = (int) (pts[idx]); + if (px < 0 || px >= w-1) + continue; + int py = (int) (pts[idx+1]); + if (py < 0 || py >= h-1) + continue; + int pz = (int) (pts[idx+2]); + if (pz < 0 || pz >= d-1) + continue; + float frac_x = pts[idx ] - px; + float frac_y = pts[idx+1] - py; + float frac_z = pts[idx+2] - pz; + + + accu[px + w*(py+h*pz)] += (1-frac_x)*(1-frac_y)*(1-frac_z) * weights[i]; + + accu[px+1 + w*(py+h*pz)] += (frac_x)*(1-frac_y)*(1-frac_z) * weights[i]; + accu[px + w*(py+1+h*pz)] += (1-frac_x)*(frac_y)*(1-frac_z) * weights[i]; + accu[px + w*(py+h*pz+h)] += (1-frac_x)*(1-frac_y)*(frac_z) * weights[i]; + + accu[px + w*(py+1+h*pz+h)] += (1-frac_x)*(frac_y)*(frac_z) * weights[i]; + accu[px+1 + w*(py+h*pz+h)] += (frac_x)*(1-frac_y)*(frac_z) * weights[i]; + accu[px+1 + w*(py+1+h*pz)] += (frac_x)*(frac_y)*(1-frac_z) * weights[i]; + + accu[px+1 + w*(py+1+h*pz+h)] += (frac_x)*(frac_y)*(frac_z) * weights[i]; + + } + + +} + + + + + + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/BuildFibres.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/BuildFibres.cpp new file mode 100644 index 0000000000..ebe5243717 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/BuildFibres.cpp @@ -0,0 +1,364 @@ +#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 +{ +public: + Particle *particles; + int pcnt; + int attrcnt; + typedef vector< Particle* > ParticleContainerType; + typedef vector< ParticleContainerType* > FiberContainerType; + + FiberContainerType* m_FiberContainer; + + CCAnalysis(float *points, int numPoints, double spacing[]) + { + m_FiberContainer = new FiberContainerType(); + particles = (Particle*) malloc(sizeof(Particle)*numPoints); + pcnt = numPoints; + attrcnt = 10; + for (int k = 0; k < numPoints; k++) + { + Particle *p = &(particles[k]); + p->R = pVector(points[attrcnt*k]/spacing[0], points[attrcnt*k+1]/spacing[1],points[attrcnt*k+2]/spacing[2]); + p->N = pVector(points[attrcnt*k+3],points[attrcnt*k+4],points[attrcnt*k+5]); + p->cap = points[attrcnt*k+6]; + p->len = points[attrcnt*k+7]; + p->mID = (int) points[attrcnt*k+8]; + p->pID = (int) points[attrcnt*k+9]; + p->ID = k; + p->label = 0; + } + } + + ~CCAnalysis() + { + for (int i=0; isize(); i++) + delete(m_FiberContainer->at(i)); + delete(m_FiberContainer); + free(particles); + } + + int 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(); + 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); + numFibers++; + } + } + } + return numFibers; + } + + void labelPredecessors(Particle *dp, ParticleContainerType* 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); + } + void labelSuccessors(Particle *dp, ParticleContainerType* container) + { + if(container->back()->ID!=dp->ID) + container->push_back(dp); + + 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/GlobalTracking/EnergyComputerBase.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/EnergyComputerBase.cpp new file mode 100644 index 0000000000..9a2f4c5711 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/EnergyComputerBase.cpp @@ -0,0 +1,368 @@ +#ifndef _ENCOMPINTERFACE +#define _ENCOMPINTERFACE + + +#include "SphereInterpolator.cpp" +#include "ParticleGrid.cpp" +#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/GlobalTracking/EnergyComputer_connec.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/EnergyComputer_connec.cpp new file mode 100644 index 0000000000..8c40d4dd43 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/EnergyComputer_connec.cpp @@ -0,0 +1,226 @@ + + + +/// bessel for wid = 1 + +//#define mbesseli0(x) ((x>=1.0) ? ((-0.2578*x+0.7236)*exp(x)) : ((x>=0.9) ? ((-0.2740*x+0.7398)*exp(x)) : ((x>=0.8) ? ((-0.3099*x+0.7720)*exp(x)) : ((x>=0.7) ? ((-0.3634*x+0.8149)*exp(x)) : ((x>=0.5) ? ((-0.4425*x+0.8663)*exp(x)) : ((x>=0.3) ? ((-0.5627*x+0.9264)*exp(x)) : ((x>=0.2) ? ((-0.6936*x+0.9657)*exp(x)) : ((x>=0.1) ? ((-0.8016*x+0.9873)*exp(x)) : ((x>=0.0) ? ((-0.9290*x+1.0000)*exp(x)) : 1 ))))))))) +//#define mbesseli0(x) ((x>=1.0) ? ((0.5652*x+0.7009)) : ((x>=0.8) ? ((0.4978*x+0.7683)) : ((x>=0.6) ? ((0.3723*x+0.8686)) : ((x>=0.4) ? ((0.2582*x+0.9371)) : ((x>=0.2) ? ((0.1519*x+0.9796)) : ((x>=0.0) ? ((0.0501*x+1.0000)) : 1 )))))) + + +inline float mbesseli0(float x) +{ + float y = x*x; + float erg = BESSEL_APPROXCOEFF[0]; + erg += y*BESSEL_APPROXCOEFF[1]; + erg += y*y*BESSEL_APPROXCOEFF[2]; + erg += y*y*y*BESSEL_APPROXCOEFF[3]; + return erg; +} + +// +// +// inline REAL mbesseli0(REAL x) +// { +// REAL y = x*x; +// REAL erg = BESSEL_APPROXCOEFF[0]; +// erg += y*BESSEL_APPROXCOEFF[1]; +// erg += y*x*BESSEL_APPROXCOEFF[2]; +// erg += x*x*BESSEL_APPROXCOEFF[3]; +// return erg; +// } + +inline float mexp(float x) +{ + + return((x>=7.0) ? 0 : ((x>=5.0) ? (-0.0029*x+0.0213) : ((x>=3.0) ? (-0.0215*x+0.1144) : ((x>=2.0) ? (-0.0855*x+0.3064) : ((x>=1.0) ? (-0.2325*x+0.6004) : ((x>=0.5) ? (-0.4773*x+0.8452) : ((x>=0.0) ? (-0.7869*x+1.0000) : 1 ))))))); + // return exp(-x); + +} + + +#include "ParticleGrid.cpp" + +#include "EnergyComputerBase.cpp" +#include + + +class EnergyComputer : public EnergyComputerBase +{ + +public: + + + float eigencon_energy; + + float chempot2; + float meanval_sq; + + float gamma_s; + float gamma_reg_s; + + float particle_weight; + float ex_strength; + float in_strength; + + float particle_length_sq; + float curv_hard; + ofstream ee_file, ie_file; + + + EnergyComputer(float *data, const int *dsz, double *cellsize, SphereInterpolator *sp, ParticleGrid *pcon, float *spimg, int spmult, vnl_matrix_fixed rotMatrix) : EnergyComputerBase(data,dsz,cellsize,sp,pcon,spimg,spmult,rotMatrix) + { +// ee_file.open("eeFile.txt"); +// ie_file.open("ieFile.txt"); + } + +// ~EnergyComputer(){ +// ee_file.close(); +// ie_file.close(); +// } + + void setParameters(float pwei,float pwid,float chempot_connection, float length,float curv_hardthres, float inex_balance, float chempot2) + { + this->chempot2 = chempot2; + + eigencon_energy = chempot_connection; + eigen_energy = 0; + particle_weight = pwei; + + float bal = 1/(1+exp(-inex_balance)); + ex_strength = 2*bal; // cleanup (todo) + in_strength = 2*(1-bal)/length/length; // cleanup (todo) + // in_strength = 0.64/length/length; // cleanup (todo) + + particle_length_sq = length*length; + curv_hard = curv_hardthres; + + float sigma_s = pwid; + gamma_s = 1/(sigma_s*sigma_s); + gamma_reg_s =1/(length*length/4); + } + + + + //////////////////////////////////////////////////////////////////////////// + ////// External Energy + //////////////////////////////////////////////////////////////////////////// + inline float computeExternalEnergy(pVector &R, pVector &N, float &cap, float &len, Particle *dp) + { + float m = SpatProb(R); + if (m == 0) + { + return -INFINITY; + } + + float Dn = evaluateODF(R,N,len); + + float Sn = 0; + float Pn = 0; + + m_ParticleGrid->computeNeighbors(R); + for (;;) + { + Particle *p = m_ParticleGrid->getNextNeighbor(); + if (p == 0) break; + if (dp != p) + { + float dot = fabs(N*p->N); + float bw = mbesseli0(dot); + float dpos = (p->R-R).norm_square(); + float w = mexp(dpos*gamma_s); + Sn += w*(bw+chempot2)*p->cap ; + w = mexp(dpos*gamma_reg_s); + Pn += w*bw; + } + } + + float energy = 0; + energy += (2*(Dn/particle_weight-Sn) - (mbesseli0(1.0)+chempot2)*cap)*cap; + + return energy*ex_strength; + } + + + //////////////////////////////////////////////////////////////////////////// + ////// Internal Energy + //////////////////////////////////////////////////////////////////////////// + + inline float computeInternalEnergy(Particle *dp) + { + float energy = eigen_energy; + + if (dp->pID != -1) + energy += computeInternalEnergyConnection(dp,+1); + + if (dp->mID != -1) + energy += computeInternalEnergyConnection(dp,-1); + + //ie_file << energy << "\n"; + + return energy; + } + + inline float computeInternalEnergyConnection(Particle *p1,int ep1) + { + Particle *p2 = 0; + int ep2; + if (ep1 == 1) + p2 = m_ParticleGrid->ID_2_address[p1->pID]; + else + p2 = m_ParticleGrid->ID_2_address[p1->mID]; + if (p2->mID == p1->ID) + ep2 = -1; + else if (p2->pID == p1->ID) + ep2 = 1; + else + fprintf(stderr,"EnergyComputer_connec: Connections are inconsistent!\n"); + + if (p2 == 0) + fprintf(stderr,"bug2"); + + return computeInternalEnergyConnection(p1,ep1,p2,ep2); + } + + inline float computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) + { +#ifdef TIMING + tic(&internalenergy_time); +#endif + + if ((p1->N*p2->N)*ep1*ep2 > -curv_hard) + return -INFINITY; + + pVector R1 = p1->R + (p1->N * (p1->len * ep1)); + pVector R2 = p2->R + (p2->N * (p2->len * ep2)); + + if ((R1-R2).norm_square() > particle_length_sq) + return -INFINITY; + + pVector R = (p2->R + p1->R)*0.5; + + if (SpatProb(R) == 0) + return -INFINITY; + + float norm1 = (R1-R).norm_square(); + float norm2 = (R2-R).norm_square(); + + + float energy = (eigencon_energy-norm1-norm2)*in_strength; + +#ifdef TIMING + toc(&internalenergy_time); +#endif + + return energy; + } + + + + + + + + + +}; + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/MersenneTwister.h b/Modules/DiffusionImaging/Tractography/GlobalTracking/MersenneTwister.h new file mode 100644 index 0000000000..d039f24113 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/MersenneTwister.h @@ -0,0 +1,594 @@ +// MersenneTwister.h +// Mersenne Twister random number generator -- a C++ class MTRand +// Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus +// Richard J. Wagner v1.0 15 May 2003 rjwagner@writeme.com + +// The Mersenne Twister is an algorithm for generating random numbers. It +// was designed with consideration of the flaws in various other generators. +// The period, 2^19937-1, and the order of equidistribution, 623 dimensions, +// are far greater. The generator is also fast; it avoids multiplication and +// division, and it benefits from caches and pipelines. For more information +// see the inventors' web page at http://www.math.keio.ac.jp/~matumoto/emt.html + +// Reference +// M. Matsumoto and T. Nishimura, "Mersenne Twister: A 623-Dimensionally +// Equidistributed Uniform Pseudo-Random Number Generator", ACM Transactions on +// Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, pp 3-30. + +// Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura, +// Copyright (C) 2000 - 2003, Richard J. Wagner +// All rights reserved. +// +// Redistribution and use in source and binary forms, with or without +// modification, are permitted provided that the following conditions +// are met: +// +// 1. Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. +// +// 2. Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in the +// documentation and/or other materials provided with the distribution. +// +// 3. The names of its contributors may not be used to endorse or promote +// products derived from this software without specific prior written +// permission. +// +// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR +// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +// The original code included the following notice: +// +// When you use this, send an email to: matumoto@math.keio.ac.jp +// with an appropriate reference to your work. +// +// It would be nice to CC: rjwagner@writeme.com and Cokus@math.washington.edu +// when you write. + +#ifndef MERSENNETWISTER_H +#define MERSENNETWISTER_H + +// Not thread safe (unless auto-initialization is avoided and each thread has +// its own MTRand object) + +#include +#include +#include +#include +#include + +class MTRand { +// Data +public: + typedef unsigned long uint32; // unsigned integer type, at least 32 bits + + enum { N = 624 }; // length of state vector + enum { SAVE = N + 1 }; // length of array for save() + +// protected: + enum { M = 397 }; // period parameter + + uint32 state[N]; // internal state + uint32 *pNext; // next value to get from state + int left; // number of values left before reload needed + + +//Methods +public: + MTRand( const uint32& oneSeed ); // initialize with a simple uint32 + MTRand( uint32 *const bigSeed, uint32 const seedLength = N ); // or an array + MTRand(); // auto-initialize with /dev/urandom or time() and clock() + ~MTRand(); + + // Do NOT use for CRYPTOGRAPHY without securely hashing several returned + // values together, otherwise the generator state can be learned after + // reading 624 consecutive values. + + // Access to 32-bit random numbers + double rand(); // real number in [0,1] + float frand(); // real number in [0,1] + double rand( const double& n ); // real number in [0,n] + double randExc(); // real number in [0,1) + double randExc( const double& n ); // real number in [0,n) + double randDblExc(); // real number in (0,1) + double randDblExc( const double& n ); // real number in (0,n) + uint32 randInt(); // integer in [0,2^32-1] + uint32 randInt( const uint32& n ); // integer in [0,n] for n < 2^32 + double operator()() { return rand(); } // same as rand() + + // Access to 53-bit random numbers (capacity of IEEE double precision) + double rand53(); // real number in [0,1) + + // Access to nonuniform random number distributions + double randNorm( const double& mean = 0.0, const double& variance = 0.0 ); + float frandn(); + + // Re-seeding functions with same behavior as initializers + void seed( const uint32 oneSeed ); + void seed( uint32 *const bigSeed, const uint32 seedLength = N ); + void seed(); + + + ////// + int Poisson(); + void initPoisson(float mean, int accuracy); + float *cumulPoisson; + int accuracyPoisson; + + int drawGamma(); + void initGamma(float theta, int shape); + float *cumulGamma; + int accuracyGamma; + + + + // Saving and loading generator state + void save( uint32* saveArray ) const; // to array of size SAVE + void load( uint32 *const loadArray ); // from such array + friend std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ); + friend std::istream& operator>>( std::istream& is, MTRand& mtrand ); + +protected: + void initialize( const uint32 oneSeed ); + void reload(); + uint32 hiBit( const uint32& u ) const { return u & 0x80000000UL; } + uint32 loBit( const uint32& u ) const { return u & 0x00000001UL; } + uint32 loBits( const uint32& u ) const { return u & 0x7fffffffUL; } + uint32 mixBits( const uint32& u, const uint32& v ) const + { return hiBit(u) | loBits(v); } + uint32 twist( const uint32& m, const uint32& s0, const uint32& s1 ) const + { return m ^ (mixBits(s0,s1)>>1) ^ (-loBit(s1) & 0x9908b0dfUL); } + static uint32 hash( time_t t, clock_t c ); +}; + + +inline MTRand::MTRand( const uint32& oneSeed ) + { seed(oneSeed); + cumulPoisson = 0; + cumulGamma = 0; + } + +inline MTRand::MTRand( uint32 *const bigSeed, const uint32 seedLength ) + { seed(bigSeed,seedLength); + cumulPoisson = 0; + cumulGamma = 0; + + } + +inline MTRand::MTRand() + { seed(); + cumulPoisson = 0; + cumulGamma = 0; + } + +inline MTRand::~MTRand() +{ + if (cumulPoisson != 0) + free(cumulPoisson); + if (cumulGamma != 0) + free(cumulGamma); +} +inline double MTRand::rand() + { return double(randInt()) * (1.0/4294967295.0); } + +inline float MTRand::frand() + { + return float(randInt()) * (1.0/4294967295.0); + } + + +inline double MTRand::rand( const double& n ) + { return rand() * n; } + +inline double MTRand::randExc() + { return double(randInt()) * (1.0/4294967296.0); } + +inline double MTRand::randExc( const double& n ) + { return randExc() * n; } + +inline double MTRand::randDblExc() + { return ( double(randInt()) + 0.5 ) * (1.0/4294967296.0); } + +inline double MTRand::randDblExc( const double& n ) + { return randDblExc() * n; } + +inline double MTRand::rand53() +{ + uint32 a = randInt() >> 5, b = randInt() >> 6; + return ( a * 67108864.0 + b ) * (1.0/9007199254740992.0); // by Isaku Wada +} + +inline double MTRand::randNorm( const double& mean, const double& variance ) +{ + // Return a real number from a normal (Gaussian) distribution with given + // mean and variance by Box-Muller method + double r = sqrt( -2.0 * log( 1.0-randDblExc()) ) * variance; + double phi = 2.0 * 3.14159265358979323846264338328 * randExc(); + return mean + r * cos(phi); +} + +inline float MTRand::frandn() +{ + // Return a real number from a normal (Gaussian) distribution with given + // mean and variance by Box-Muller method + float r = sqrt( -2.0 * log( 1.0-randDblExc()) ); + double phi = 2.0 * 3.14159265358979323846264338328 * randExc(); + return r * cos(phi); +} + +inline int MTRand::Poisson() +{ + int j; + + float r = frand(); + + int rl = 0; + int rh = accuracyPoisson-1; + while(rh != rl) + { + j = rl + (rh-rl)/2; + if (r < cumulPoisson[j]) + { + rh = j; + continue; + } + if (r > cumulPoisson[j]) + { + rl = j+1; + continue; + } + break; + } + j = rh; + + return j; +} + + +inline int MTRand::drawGamma() +{ + int j; + + float r = frand(); + + int rl = 0; + int rh = accuracyGamma-1; + while(rh != rl) + { + j = rl + (rh-rl)/2; + if (r < cumulGamma[j]) + { + rh = j; + continue; + } + if (r > cumulGamma[j]) + { + rl = j+1; + continue; + } + break; + } + j = rh; + + return j; +} + + +/* +inline void MTRand::initPoisson(float mean,int accuracy) +{ + float p[accuracy]; + float Z = exp(-mean); + p[0] = 1; + for (int i = 1; i < accuracy;i++) + p[i] = p[i-1]*mean/i; + for (int i = 0; i < accuracy;i++) + p[i] *= Z; + + if (cumulPoisson != 0) + free(cumulPoisson); + + cumulPoisson = (float*) malloc(sizeof(float)*accuracy); + cumulPoisson[0] = p[0]; + for (int i = 1; i < accuracy;i++) + cumulPoisson[i] = cumulPoisson[i-1]+p[i]; + for (int i = 0; i < accuracy;i++) + cumulPoisson[i] = cumulPoisson[i]/cumulPoisson[accuracy-1]; + + accuracyPoisson = accuracy; + + +} + + + + + +inline void MTRand::initGamma(float theta,int k) +{ + int accuracy = int(k*theta*5); + + float p[accuracy]; + int fac = 1; + for (int i = 1; i < k-1;i++) + fac *= i; + + for (int i = 0; i < accuracy;i++) + { + p[i] = pow(i/theta,k-1)/fac * exp(i/theta); + } + if (cumulGamma != 0) + free(cumulGamma); + + cumulGamma = (float*) malloc(sizeof(float)*accuracy); + cumulGamma[0] = p[0]; + for (int i = 1; i < accuracy;i++) + cumulGamma[i] = cumulGamma[i-1]+p[i]; + for (int i = 0; i < accuracy;i++) + cumulGamma[i] = cumulGamma[i]/cumulGamma[accuracy-1]; + + accuracyGamma = accuracy; + + +} + + +*/ + + + + +inline MTRand::uint32 MTRand::randInt() +{ + // Pull a 32-bit integer from the generator state + // Every other access function simply transforms the numbers extracted here + + if( left == 0 ) reload(); + --left; + + register uint32 s1; + s1 = *pNext++; + s1 ^= (s1 >> 11); + s1 ^= (s1 << 7) & 0x9d2c5680UL; + s1 ^= (s1 << 15) & 0xefc60000UL; + return ( s1 ^ (s1 >> 18) ); +} + +inline MTRand::uint32 MTRand::randInt( const uint32& n ) +{ + // Find which bits are used in n + // Optimized by Magnus Jonsson (magnus@smartelectronix.com) + uint32 used = n; + used |= used >> 1; + used |= used >> 2; + used |= used >> 4; + used |= used >> 8; + used |= used >> 16; + + // Draw numbers until one is found in [0,n] + uint32 i; + do + i = randInt() & used; // toss unused bits to shorten search + while( i > n ); + return i; +} + + +inline void MTRand::seed( const uint32 oneSeed ) +{ + // Seed the generator with a simple uint32 + initialize(oneSeed); + reload(); +} + + +inline void MTRand::seed( uint32 *const bigSeed, const uint32 seedLength ) +{ + // Seed the generator with an array of uint32's + // There are 2^19937-1 possible initial states. This function allows + // all of those to be accessed by providing at least 19937 bits (with a + // default seed length of N = 624 uint32's). Any bits above the lower 32 + // in each element are discarded. + // Just call seed() if you want to get array from /dev/urandom + initialize(19650218UL); + register int i = 1; + register uint32 j = 0; + register int k = ( N > seedLength ? N : seedLength ); + for( ; k; --k ) + { + state[i] = + state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1664525UL ); + state[i] += ( bigSeed[j] & 0xffffffffUL ) + j; + state[i] &= 0xffffffffUL; + ++i; ++j; + if( i >= N ) { state[0] = state[N-1]; i = 1; } + if( j >= seedLength ) j = 0; + } + for( k = N - 1; k; --k ) + { + state[i] = + state[i] ^ ( (state[i-1] ^ (state[i-1] >> 30)) * 1566083941UL ); + state[i] -= i; + state[i] &= 0xffffffffUL; + ++i; + if( i >= N ) { state[0] = state[N-1]; i = 1; } + } + state[0] = 0x80000000UL; // MSB is 1, assuring non-zero initial array + reload(); +} + + +inline void MTRand::seed() +{ + // Seed the generator with an array from /dev/urandom if available + // Otherwise use a hash of time() and clock() values + + // First try getting an array from /dev/urandom + FILE* urandom = fopen( "/dev/urandom", "rb" ); + if( urandom ) + { + uint32 bigSeed[N]; + register uint32 *s = bigSeed; + register int i = N; + register bool success = true; + while( success && i-- ) + success = fread( s++, sizeof(uint32), 1, urandom ); + fclose(urandom); + if( success ) { seed( bigSeed, N ); return; } + } + + // Was not successful, so use time() and clock() instead + seed( hash( time(NULL), clock() ) ); +} + + +inline void MTRand::initialize( const uint32 seed ) +{ + // Initialize generator state with seed + // See Knuth TAOCP Vol 2, 3rd Ed, p.106 for multiplier. + // In previous versions, most significant bits (MSBs) of the seed affect + // only MSBs of the state array. Modified 9 Jan 2002 by Makoto Matsumoto. + register uint32 *s = state; + register uint32 *r = state; + register int i = 1; + *s++ = seed & 0xffffffffUL; + for( ; i < N; ++i ) + { + *s++ = ( 1812433253UL * ( *r ^ (*r >> 30) ) + i ) & 0xffffffffUL; + r++; + } +} + + +inline void MTRand::reload() +{ + // Generate N new values in state + // Made clearer and faster by Matthew Bellew (matthew.bellew@home.com) + register uint32 *p = state; + register int i; + for( i = N - M; i--; ++p ) + *p = twist( p[M], p[0], p[1] ); + for( i = M; --i; ++p ) + *p = twist( p[M-N], p[0], p[1] ); + *p = twist( p[M-N], p[0], state[0] ); + + left = N, pNext = state; +} + + +inline MTRand::uint32 MTRand::hash( time_t t, clock_t c ) +{ + // Get a uint32 from t and c + // Better than uint32(x) in case x is floating point in [0,1] + // Based on code by Lawrence Kirby (fred@genesis.demon.co.uk) + + static uint32 differ = 0; // guarantee time-based seeds will change + + uint32 h1 = 0; + unsigned char *p = (unsigned char *) &t; + for( size_t i = 0; i < sizeof(t); ++i ) + { + h1 *= UCHAR_MAX + 2U; + h1 += p[i]; + } + uint32 h2 = 0; + p = (unsigned char *) &c; + for( size_t j = 0; j < sizeof(c); ++j ) + { + h2 *= UCHAR_MAX + 2U; + h2 += p[j]; + } + return ( h1 + differ++ ) ^ h2; +} + + +inline void MTRand::save( uint32* saveArray ) const +{ + register uint32 *sa = saveArray; + register const uint32 *s = state; + register int i = N; + for( ; i--; *sa++ = *s++ ) {} + *sa = left; +} + + +inline void MTRand::load( uint32 *const loadArray ) +{ + register uint32 *s = state; + register uint32 *la = loadArray; + register int i = N; + for( ; i--; *s++ = *la++ ) {} + left = *la; + pNext = &state[N-left]; +} + + +inline std::ostream& operator<<( std::ostream& os, const MTRand& mtrand ) +{ + register const MTRand::uint32 *s = mtrand.state; + register int i = mtrand.N; + for( ; i--; os << *s++ << "\t" ) {} + return os << mtrand.left; +} + + +inline std::istream& operator>>( std::istream& is, MTRand& mtrand ) +{ + register MTRand::uint32 *s = mtrand.state; + register int i = mtrand.N; + for( ; i--; is >> *s++ ) {} + is >> mtrand.left; + mtrand.pNext = &mtrand.state[mtrand.N-mtrand.left]; + return is; +} + +#endif // MERSENNETWISTER_H + +// Change log: +// +// v0.1 - First release on 15 May 2000 +// - Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus +// - Translated from C to C++ +// - Made completely ANSI compliant +// - Designed convenient interface for initialization, seeding, and +// obtaining numbers in default or user-defined ranges +// - Added automatic seeding from /dev/urandom or time() and clock() +// - Provided functions for saving and loading generator state +// +// v0.2 - Fixed bug which reloaded generator one step too late +// +// v0.3 - Switched to clearer, faster reload() code from Matthew Bellew +// +// v0.4 - Removed trailing newline in saved generator format to be consistent +// with output format of built-in types +// +// v0.5 - Improved portability by replacing static const int's with enum's and +// clarifying return values in seed(); suggested by Eric Heimburg +// - Removed MAXINT constant; use 0xffffffffUL instead +// +// v0.6 - Eliminated seed overflow when uint32 is larger than 32 bits +// - Changed integer [0,n] generator to give better uniformity +// +// v0.7 - Fixed operator precedence ambiguity in reload() +// - Added access for real numbers in (0,1) and (0,n) +// +// v0.8 - Included time.h header to properly support time_t and clock_t +// +// v1.0 - Revised seeding to match 26 Jan 2002 update of Nishimura and Matsumoto +// - Allowed for seeding with arrays of any length +// - Added access for real numbers in [0,1) with 53-bit resolution +// - Added access for real numbers from normal (Gaussian) distributions +// - Increased overall speed by optimizing twist() +// - Doubled speed of integer [0,n] generation +// - Fixed out-of-range number generation on 64-bit machines +// - Improved portability by substituting literal constants for long enum's +// - Changed license from GNU LGPL to BSD diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/ParticleGrid.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/ParticleGrid.cpp new file mode 100644 index 0000000000..babd7b47fd --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/ParticleGrid.cpp @@ -0,0 +1,611 @@ + + + +#ifndef _PARTICLEGRID +#define _PARTICLEGRID + + +#include "auxilary_classes.cpp" + + +template + class ParticleGrid +{ + + //////////////// Container +public: + T *particles; // particles in linear array + int pcnt; // actual number of particles + int concnt; // number of connections + int celloverflows; + + T **ID_2_address; + +private: + + int capacity; // maximal number of particles + int increase_step; + + /////////////////// Grid + + T **grid; // the grid + + // grid size + int nx; + int ny; + int nz; + + // scaling factor for grid + float mulx; + float muly; + float mulz; + + int csize; // particle capacity of single cell in grid + int *occnt; // occupation count of grid cells + int gridsize; // total number of cells + float m_Memory; + + struct NeighborTracker // to run over the neighbors + { + int cellidx[8]; + int cellidx_c[8]; + + int cellcnt; + int pcnt; + + } nbtrack; + + + +public: + + + ParticleGrid() + { + + //// involving the container + capacity = 0; + particles = 0; + ID_2_address = 0; + pcnt = 0; + concnt = 0; + celloverflows = 0; + + ////// involvin the grid + nx = 0; ny = 0; nz = 0; csize = 0; + gridsize = 0; + grid = (T**) 0; + occnt = (int*) 0; + + increase_step = 100000; + m_Memory = 0; + } + + float GetMemoryUsage() + { + return m_Memory; + } + + int allocate(int _capacity, + int _nx, int _ny, int _nz, float cellsize, int cellcapacity) + { + //// involving the container + capacity = _capacity; + particles = (T*) malloc(sizeof(T)*capacity); + ID_2_address = (T**) malloc(sizeof(T*)*capacity); + + if (particles == 0 || ID_2_address == 0) + { + fprintf(stderr,"error: Out of Memory\n"); + capacity = 0; + return -1; + } + else + { + fprintf(stderr,"Allocated Memory for %i particles \n",capacity); + } + + pcnt = 0; + int i; + for (i = 0;i < capacity;i++) + { + ID_2_address[i] = &(particles[i]); // initialize pointer in LUT + particles[i].ID = i; // initialize unique IDs + } + + ////// involvin the grid + nx = _nx; ny = _ny; nz = _nz; csize = cellcapacity; + gridsize = nx*ny*nz; + + m_Memory = (float)(sizeof(T*)*gridsize*csize)/1000000; + + grid = (T**) malloc(sizeof(T*)*gridsize*csize); + occnt = (int*) malloc(sizeof(int)*gridsize); + + if (grid == 0 || occnt == 0) + { + fprintf(stderr,"error: Out of Memory\n"); + capacity = 0; + return -1; + } + + for (i = 0;i < gridsize;i++) + occnt[i] = 0; + + mulx = 1/cellsize; + muly = 1/cellsize; + mulz = 1/cellsize; + + return 1; + } + + + + int reallocate() + { + int new_capacity = capacity + increase_step; + T* new_particles = (T*) realloc(particles,sizeof(T)*new_capacity); + T** new_ID_2_address = (T**) realloc(ID_2_address,sizeof(T*)*new_capacity); + if (new_particles == 0 || new_ID_2_address == 0) + { + fprintf(stderr,"ParticleGird:reallocate: out of memory!\n"); + return -1; + } + fprintf(stderr," now %i particles are allocated \n",new_capacity); + m_Memory = (float)(sizeof(T*)*new_capacity)/1000000; + + int i; + for (i = 0; i < capacity; i++) + { + new_ID_2_address[i] = new_ID_2_address[i] - particles + new_particles; // shift address + } + for (i = capacity; i < new_capacity; i++) + { + new_particles[i].ID = i; // initialize unique IDs + new_ID_2_address[i] = &(new_particles[i]) ; // initliaze pointer in LUT + } + for (i = 0; i < nx*ny*nz*csize; i++) + { + grid[i] = grid[i] - particles + new_particles; + } + particles = new_particles; + ID_2_address = new_ID_2_address; + capacity = new_capacity; + + return 1; + } + + + + ~ParticleGrid() + { + if (particles != 0) + free(particles); + if (grid != 0) + free(grid); + if (occnt != 0) + free(occnt); + if (ID_2_address != 0) + free(ID_2_address); + } + + + + int ID_2_index(int ID) + { + if (ID == -1) + return -1; + else + return (ID_2_address[ID] - particles); + + } + + + T* newParticle(pVector R) + { + /////// get free place in container; + if (pcnt >= capacity) + { + fprintf(stderr,"capacity overflow , reallocating ...\n"); + if (reallocate() == -1) + { + fprintf(stderr,"out of Memory!!\n"); + return 0; + } + } + + int xint = int(R[0]*mulx); + if (xint < 0) { //fprintf(stderr,"error: out of grid\n"); + return 0;} + if (xint >= nx) { // fprintf(stderr,"error: out of grid\n"); + return 0;} + int yint = int(R[1]*muly); + if (yint < 0) { //fprintf(stderr,"error: out of grid\n"); + return 0;} + if (yint >= ny) {// fprintf(stderr,"error: out of grid\n"); + return 0;} + int zint = int(R[2]*mulz); + if (zint < 0) {// fprintf(stderr,"error: out of grid\n"); + return 0;} + if (zint >= nz) { //fprintf(stderr,"error: out of grid\n"); + return 0;} + + int idx = xint + nx*(yint + ny*zint); + if (occnt[idx] < csize) + { + T *p = &(particles[pcnt]); + p->R = R; + p->mID = -1; + p->pID = -1; + pcnt++; + p->gridindex = csize*idx + occnt[idx]; + grid[p->gridindex] = p; + occnt[idx]++; + return p; + } + else + { + celloverflows++; + //fprintf(stderr,"error: cell overflow \n"); + return 0; + } + } + + + inline void updateGrid(int k) + { + T* p = &(particles[k]); + + /////// find new grid cell + int xint = int(p->R[0]*mulx); + if (xint < 0) { remove(k); return; } + if (xint >= nx) { remove(k); return; } + int yint = int(p->R[1]*muly); + if (yint < 0) { remove(k); return; } + if (yint >= ny) { remove(k); return; } + int zint = int(p->R[2]*mulz); + if (zint < 0) { remove(k); return; } + if (zint >= nz) { remove(k); return; } + + + int idx = xint + nx*(yint+ zint*ny); + int cellidx = p->gridindex/csize; + if (idx != cellidx) // cell has changed + { + + if (occnt[idx] < csize) + { + // remove from old position in grid; + int grdindex = p->gridindex; + grid[grdindex] = grid[cellidx*csize + occnt[cellidx]-1]; + grid[grdindex]->gridindex = grdindex; + occnt[cellidx]--; + + // insert at new position in grid + p->gridindex = idx*csize + occnt[idx]; + grid[p->gridindex] = p; + occnt[idx]++; + } + else + remove(k); + + } + } + + + inline bool tryUpdateGrid(int k) + { + T* p = &(particles[k]); + + /////// find new grid cell + int xint = int(p->R[0]*mulx); + if (xint < 0) { return false; } + if (xint >= nx) { return false; } + int yint = int(p->R[1]*muly); + if (yint < 0) { return false; } + if (yint >= ny) { return false; } + int zint = int(p->R[2]*mulz); + if (zint < 0) { return false; } + if (zint >= nz) { return false; } + + + int idx = xint + nx*(yint+ zint*ny); + int cellidx = p->gridindex/csize; + if (idx != cellidx) // cell has changed + { + + if (occnt[idx] < csize) + { + // remove from old position in grid; + int grdindex = p->gridindex; + grid[grdindex] = grid[cellidx*csize + occnt[cellidx]-1]; + grid[grdindex]->gridindex = grdindex; + occnt[cellidx]--; + + // insert at new position in grid + p->gridindex = idx*csize + occnt[idx]; + grid[p->gridindex] = p; + occnt[idx]++; + return true; + } + else + return false; + + } + return true; + } + + + + inline void remove(int k) + { + T* p = &(particles[k]); + int grdindex = p->gridindex; + int cellidx = grdindex/csize; + int idx = grdindex%csize; + + // remove pending connections + if (p->mID != -1) + destroyConnection(p,-1); + if (p->pID != -1) + destroyConnection(p,+1); + + // remove from grid + if (idx < occnt[cellidx]-1) + { + grid[grdindex] = grid[cellidx*csize + occnt[cellidx]-1]; + grid[grdindex]->gridindex = grdindex; + } + occnt[cellidx]--; + + + + // remove from container + if (kID; + int move_ID = particles[pcnt-1].ID; + + *p = particles[pcnt-1]; // move very last particle to empty slot + particles[pcnt-1].ID = todel_ID; // keep IDs unique + grid[p->gridindex] = p; // keep gridindex consistent + + // permute address table + ID_2_address[todel_ID] = &(particles[pcnt-1]); + ID_2_address[move_ID] = p; + + } + pcnt--; + + + + } + + inline void computeNeighbors(pVector &R) + { + float xfrac = R.GetX()*mulx; + float yfrac = R.GetY()*muly; + float zfrac = R.GetZ()*mulz; + int xint = int(xfrac); + int yint = int(yfrac); + int zint = int(zfrac); + + int dx = -1; + if (xfrac-xint > 0.5) dx = 1; + if (xint <= 0) { xint = 0; dx = 1; } + if (xint >= nx-1) { xint = nx-1; dx = -1; } + + int dy = -1; + if (yfrac-yint > 0.5) dy = 1; + if (yint <= 0) {yint = 0; dy = 1; } + if (yint >= ny-1) {yint = ny-1; dy = -1;} + + int dz = -1; + if (zfrac-zint > 0.5) dz = 1; + if (zint <= 0) {zint = 0; dz = 1; } + if (zint >= nz-1) {zint = nz-1; dz = -1;} + + + nbtrack.cellidx[0] = xint + nx*(yint+zint*ny); + nbtrack.cellidx[1] = nbtrack.cellidx[0] + dx; + nbtrack.cellidx[2] = nbtrack.cellidx[1] + dy*nx; + nbtrack.cellidx[3] = nbtrack.cellidx[2] - dx; + nbtrack.cellidx[4] = nbtrack.cellidx[0] + dz*nx*ny; + nbtrack.cellidx[5] = nbtrack.cellidx[4] + dx; + nbtrack.cellidx[6] = nbtrack.cellidx[5] + dy*nx; + nbtrack.cellidx[7] = nbtrack.cellidx[6] - dx; + + + nbtrack.cellidx_c[0] = csize*nbtrack.cellidx[0]; + nbtrack.cellidx_c[1] = csize*nbtrack.cellidx[1]; + nbtrack.cellidx_c[2] = csize*nbtrack.cellidx[2]; + nbtrack.cellidx_c[3] = csize*nbtrack.cellidx[3]; + nbtrack.cellidx_c[4] = csize*nbtrack.cellidx[4]; + nbtrack.cellidx_c[5] = csize*nbtrack.cellidx[5]; + nbtrack.cellidx_c[6] = csize*nbtrack.cellidx[6]; + nbtrack.cellidx_c[7] = csize*nbtrack.cellidx[7]; + + + + nbtrack.cellcnt = 0; + nbtrack.pcnt = 0; + + } + + + + inline T *getNextNeighbor() + { + + if (nbtrack.pcnt < occnt[nbtrack.cellidx[nbtrack.cellcnt]]) + { + + return grid[nbtrack.cellidx_c[nbtrack.cellcnt] + (nbtrack.pcnt++)]; + + } + else + { + + for(;;) + { + nbtrack.cellcnt++; + if (nbtrack.cellcnt >= 8) + return 0; + if (occnt[nbtrack.cellidx[nbtrack.cellcnt]] > 0) + break; + } + + nbtrack.pcnt = 1; + return grid[nbtrack.cellidx_c[nbtrack.cellcnt]]; + } + } + + + inline void createConnection(T *P1,int ep1, T *P2, int ep2) + { + if (ep1 == -1) + P1->mID = P2->ID; + else + P1->pID = P2->ID; + + if (ep2 == -1) + P2->mID = P1->ID; + else + P2->pID = P1->ID; + + concnt++; + } + + inline void destroyConnection(T *P1,int ep1, T *P2, int ep2) + { + if (ep1 == -1) + P1->mID = -1; + else + P1->pID = -1; + + if (ep2 == -1) + P2->mID = -1; + else + P2->pID = -1; + concnt--; + } + + inline void destroyConnection(T *P1,int ep1) + { + + T *P2 = 0; + int ep2; + if (ep1 == 1) + { + P2 = ID_2_address[P1->pID]; + P1->pID = -1; + } + else + { + P2 = ID_2_address[P1->mID]; + P1->mID = -1; + } + if (P2->mID == P1->ID) + { + P2->mID = -1; + } + else + { + P2->pID = -1; + } + concnt--; + + } +}; + + + + +/* + + +struct Connection +{ + int lID; + int rID; +}; + + +class ConnectionContainer +{ + + //////////////// Container +public: + Connection *cons; // cons in linear array + int ccnt; // actual number of cons + +private: + + int capacity; // maximal number of particles + + + +public: + + + ConnectionContainer() + { + + //// involving the container + capacity = 0; + cons = 0; + ccnt = 0; + + } + + + + void allocate(int _capacity) + { + + //// involving the container + capacity = _capacity; + cons = (Connection*) malloc(sizeof(Connection)*capacity); + ccnt = 0; + + } + + ~ConnectionContainer() + { + if (cons != 0) + free(cons); + } + + + Connection* newConnection(int lid,int rid) + { + /////// get free place in container; + if (ccnt < capacity) + { + Connection *c = &(cons[ccnt]); + c->lID = lid; + c->rID = rid; + ccnt++; + return c; + } + return 0; + } + + inline void remove(int k) + { + Connection* c = &(cons[k]); + + // remove from container + if (k +#include "MersenneTwister.h" + +class RJMCMCBase +{ +public: + + ParticleGrid m_ParticleGrid; + float *m_QBallImageData; + const int *datasz; + EnergyComputerBase *enc; + int m_Iterations; + float width; + float height; + float depth; + double *voxsize; + int m_NumAttributes; + int m_AcceptedProposals; + + RJMCMCBase(float *points,int numPoints, float *dimg, const int *dsz, double *voxsize, double cellsize) + : m_QBallImageData(dimg) + , datasz(dsz) + , enc(0) + , width(dsz[0]*voxsize[0]) + , height(dsz[1]*voxsize[1]) + , depth(dsz[2]*voxsize[2]) + , voxsize(voxsize) + , m_NumAttributes(0) + , m_AcceptedProposals(0) + { + fprintf(stderr,"Data dimensions (mm) : %f x %f x %f\n",width,height,depth); + fprintf(stderr,"Data dimensions (voxel) : %i x %i x %i\n",datasz[0],datasz[1],datasz[2]); + fprintf(stderr,"voxel size (mm) : %lf x %lf x %lf\n",voxsize[0],voxsize[1],voxsize[2]); + + float cellcnt_x = (int)((float)width/cellsize) +1; + float cellcnt_y = (int)((float)height/cellsize) +1; + float cellcnt_z = (int)((float)depth/cellsize) +1; + //int cell_capacity = 2048; + int cell_capacity = 64; + + fprintf(stderr,"grid dimensions : %f x %f x %f\n",cellcnt_x,cellcnt_y,cellcnt_z); + fprintf(stderr,"grid cell size (mm) : %f^3\n",cellsize); + fprintf(stderr,"cell capacity : %i\n",cell_capacity); + fprintf(stderr,"#cells*cellcap : %.1f K\n",cell_capacity*cellcnt_x*cellcnt_y*cellcnt_z/1000); + + int minsize = 1000000; + int err = m_ParticleGrid.allocate(((numPoints>minsize)? (numPoints+100000) : minsize), cellcnt_x, cellcnt_y, cellcnt_z, cellsize, cell_capacity); + + if (err == -1) + { + fprintf(stderr,"RJMCMCBase: out of Memory!\n"); + return; + } + + m_NumAttributes = 10; + for (int k = 0; k < numPoints; k++) + { + Particle *p = m_ParticleGrid.newParticle(pVector(points[m_NumAttributes*k], points[m_NumAttributes*k+1],points[m_NumAttributes*k+2])); + if (p!=0) + { + p->N = pVector(points[m_NumAttributes*k+3],points[m_NumAttributes*k+4],points[m_NumAttributes*k+5]); + p->cap = points[m_NumAttributes*k+6]; + p->len = points[m_NumAttributes*k+7]; + p->mID = (int) points[m_NumAttributes*k+8]; + p->pID = (int) points[m_NumAttributes*k+9]; + if (p->mID != -1) + m_ParticleGrid.concnt++; + if (p->pID != -1) + m_ParticleGrid.concnt++; + p->label = 0; + } + else + { + fprintf(stderr,"error: cannot allocate particle, con. indices will be wrong! \n"); + } + } + m_ParticleGrid.concnt /= 2; + + m_Iterations = 0; + m_AcceptedProposals = 0; + } + + ~RJMCMCBase() + { + + } + + void WriteOutParticles(float *npoints) + { + for (int k = 0; k < m_ParticleGrid.pcnt; k++) + { + Particle *p = &(m_ParticleGrid.particles[k]); + npoints[m_NumAttributes*k] = p->R.GetX(); + npoints[m_NumAttributes*k+1] = p->R.GetY(); + npoints[m_NumAttributes*k+2] = p->R.GetZ(); + npoints[m_NumAttributes*k+3] = p->N.GetX(); + npoints[m_NumAttributes*k+4] = p->N.GetY(); + npoints[m_NumAttributes*k+5] = p->N.GetZ(); + npoints[m_NumAttributes*k+6] = p->cap; + npoints[m_NumAttributes*k+7] = p->len; + npoints[m_NumAttributes*k+8] = m_ParticleGrid.ID_2_index(p->mID); + npoints[m_NumAttributes*k+9] = m_ParticleGrid.ID_2_index(p->pID); + } + } + + void SetEnergyComputer(EnergyComputerBase *e) + { + enc = e; + } + + void Iterate(float* acceptance, unsigned long* numCon, unsigned long* numPart, bool *abort) + { + m_AcceptedProposals = 0; + for (int it = 0; it < m_Iterations;it++) + { + if (*abort) + break; + + IterateOneStep(); + + *numCon = m_ParticleGrid.concnt; + *numPart = m_ParticleGrid.pcnt; + } + *acceptance = (float)m_AcceptedProposals/m_Iterations; + } + + virtual void IterateOneStep() + { + + } +}; + + + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/RJMCMC_randshift.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/RJMCMC_randshift.cpp new file mode 100644 index 0000000000..0b83f23287 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/RJMCMC_randshift.cpp @@ -0,0 +1,618 @@ + +#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); + //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/GlobalTracking/SphereInterpolator.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/SphereInterpolator.cpp new file mode 100644 index 0000000000..0f37dd624d --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/SphereInterpolator.cpp @@ -0,0 +1,122 @@ + + +#include "auxilary_classes.cpp" + + +class SphereInterpolator +{ +public: + float *barycoords; + int *indices; + int size; // size of LUT + int sN; // (sizeofLUT-1)/2 + int nverts; // number of data vertices + + + float beta; + float inva; + float b; + + + int *idx; + float *interpw; + + SphereInterpolator(float *barycoords, int *indices, int numverts, int sizeLUT, float beta) + { + this->barycoords = barycoords; + this->indices = indices; + this->size = sizeLUT; + this->sN = (sizeLUT-1)/2; + this->nverts = numverts; + this->beta = beta; + + inva = (sqrt(1+beta)-sqrt(beta)); + b = 1/(1-sqrt(1/beta + 1)); + + } + + + inline void getInterpolation(vnl_vector_fixed N) + { + float nx = N[0]; + float ny = N[1]; + float nz = N[2]; + + if (nz > 0.5) + { + int x = float2int(nx); + int y = float2int(ny); + int i = 3*6*(x+y*size); // (:,1,x,y) + idx = indices+i; + interpw = barycoords +i; + return; + } + if (nz < -0.5) + { + int x = float2int(nx); + int y = float2int(ny); + int i = 3*(1+6*(x+y*size)); // (:,2,x,y) + idx = indices+i; + interpw = barycoords +i; + return; + } + if (nx > 0.5) + { + int z = float2int(nz); + int y = float2int(ny); + int i = 3*(2+6*(z+y*size)); // (:,2,x,y) + idx = indices+i; + interpw = barycoords +i; + return; + } + if (nx < -0.5) + { + int z = float2int(nz); + int y = float2int(ny); + int i = 3*(3+6*(z+y*size)); // (:,2,x,y) + idx = indices+i; + interpw = barycoords +i; + return; + } + if (ny > 0) + { + int x = float2int(nx); + int z = float2int(nz); + int i = 3*(4+6*(x+z*size)); // (:,1,x,y) + idx = indices+i; + interpw = barycoords +i; + return; + } + else + { + int x = float2int(nx); + int z = float2int(nz); + int i = 3*(5+6*(x+z*size)); // (:,1,x,y) + idx = indices+i; + interpw = barycoords +i; + return; + } + + } + + + inline float invrescale(float f) + { + float x = (fabs(f)-b)*inva; + if (f>0) + return (x*x-beta); + else + return beta - x*x; + } + + inline int float2int(float x) + { + return int((invrescale(x)+1)*sN-0.5); + + } + + +}; + + + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/anisoDiffusion.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/anisoDiffusion.cpp new file mode 100644 index 0000000000..93e593ab51 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/anisoDiffusion.cpp @@ -0,0 +1,134 @@ + +#include "mex.h" +#include "matrix.h" +#include +#include +#define REAL float + +void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) +{ + + if(nrhs != 3) { + printf("\nUsage: res = aniosDiffusion(img,Dcoeff,[alpha it])\n\n",nrhs); + printf(" Propagates Heat equation with position dependedd conductivity Dcoeff\n\n"); + printf(" img - inital state\n"); + printf(" Dcoeff - contains Dtensor of size [ 6 size(img)] with compnents xx yy zz xz yz xz \n"); + printf(" alpha - timestep \n"); + printf(" it - number of iterations \n"); + + return; + } else if(nlhs>1) { + printf("Too many output arguments\n"); + return; + } + + int pcnt = 0; + const mxArray *Img; + Img = prhs[pcnt++]; + const int numdim = mxGetNumberOfDimensions(Img); + const int *dims = mxGetDimensions(Img); + REAL *img = (REAL*) mxGetData(Img); + if (numdim != 3) + { + mexPrintf("input has to be 3dim\n"); + return; + } + int w = dims[0]; + int h = dims[1]; + int d = dims[2]; + + const mxArray *BG; + BG = prhs[pcnt++]; + const int *dimsbg = mxGetDimensions(BG); + if (dimsbg[0] != 6) + { + mexPrintf("blurrguide has to be rank 2\n"); + return; + } + REAL *blurrg = (REAL*) mxGetData(BG); + + + const mxArray *Factors; + Factors = prhs[pcnt++]; + REAL *factors = (REAL*) mxGetData(Factors); + REAL alpha = factors[0]; + int itmax = (int) factors[1]; + + plhs[0] = mxCreateNumericArray(3,dims,mxGetClassID(Img),mxREAL); + REAL *result = (REAL*) mxGetData(plhs[0]); + memcpy(result,img,sizeof(REAL)*w*h*d); + + + + int ndims[4]; + ndims[0] = 3; ndims[1] = dims[0]; ndims[2] = dims[1]; ndims[3] = dims[2]; + const mxArray *Tmp1 = mxCreateNumericArray(4,ndims,mxGetClassID(Img),mxREAL); + REAL *tmp1 = (REAL*) mxGetData(Tmp1); + + + + int wh = w*h; + + + for (int k = 0; k < itmax ; k++) + { + + for (int z = 1; z < d-1; z++) + for (int y = 1; y < h-1; y++) + for (int x = 1; x < w-1; x++) + { + int idx = z*wh + y*w + x; + REAL gx = result[idx-1] - result[idx]; + REAL gy = result[idx-w] - result[idx]; + REAL gz = result[idx-wh] - result[idx]; + REAL *D = &(blurrg[6*idx]); + REAL *tmp = &(tmp1[3*idx]); + tmp[0] = gx*D[0] + gy*D[3] + gz*D[4]; + tmp[1] = gx*D[3] + gy*D[1] + gz*D[5]; + tmp[2] = gx*D[4] + gy*D[5] + gz*D[2]; + } + + for (int z = 1; z < d-1; z++) + for (int y = 1; y < h-1; y++) + for (int x = 1; x < w-1; x++) + { + int idx = (z*wh + y*w + x); + REAL *tmp = &(tmp1[3*idx]); + result[idx] = result[idx] - alpha*(tmp[3]-tmp[0] + tmp[1+3*w]-tmp[1] + tmp[2+3*wh]-tmp[2]); + } + /* + + + for (int z = 1; z < d-1; z++) + for (int y = 1; y < h-1; y++) + for (int x = 1; x < w-1; x++) + { + int idx = z*wh + y*w + x; + REAL gxx = result[idx+1] + result[idx-1] - 2*result[idx]; + REAL gyy = result[idx+w] + result[idx-w] - 2*result[idx]; + REAL gzz = result[idx+wh] + result[idx-wh] - 2*result[idx]; + REAL gxy = (result[idx+1+w] - result[idx-1+w] - (result[idx+1-w] - result[idx-1-w]))/4; + REAL gxz = (result[idx+1+wh] - result[idx-1+wh] - (result[idx+1-wh] - result[idx-1-wh]))/4; + REAL gyz = (result[idx+w+wh] - result[idx-w+wh] - (result[idx+w-wh] - result[idx-w-wh]))/4; + + REAL *D = &(blurrg[6*idx]); + //tmp[0] = gxx*D[0] + gxy*D[3] + gxz*D[4]; + //tmp[1] = gxy*D[3] + gyy*D[1] + gyz*D[5]; + //tmp[2] = gxz*D[4] + gyz*D[5] + gzz*D[2]; + + result[idx] = result[idx] - alpha*( gxx*D[0] + gxy*D[3] + gxz*D[4] + gxy*D[3] + gyy*D[1] + gyz*D[5] + gxz*D[4] + gyz*D[5] + gzz*D[2]); + } + */ + } + + + + + + +} + + + + + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/anisoDiffusionHomogenous.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/anisoDiffusionHomogenous.cpp new file mode 100644 index 0000000000..d5a666b138 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/anisoDiffusionHomogenous.cpp @@ -0,0 +1,134 @@ + +#include "mex.h" +#include "matrix.h" +#include +#include +#define REAL float + +void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) +{ + + if(nrhs != 3) { + printf("\nUsage: res = aniosDiffusion(img,Dcoeff,[alpha it])\n\n",nrhs); + printf(" Propagates Heat equation with position dependedd conductivity Dcoeff\n\n"); + printf(" img - inital state\n"); + printf(" Dcoeff - contains Dtensor of size [ 6 size(img)] with compnents xx yy zz xz yz xz \n"); + printf(" alpha - timestep \n"); + printf(" it - number of iterations \n"); + + return; + } else if(nlhs>1) { + printf("Too many output arguments\n"); + return; + } + + int pcnt = 0; + const mxArray *Img; + Img = prhs[pcnt++]; + const int numdim = mxGetNumberOfDimensions(Img); + const int *dims = mxGetDimensions(Img); + REAL *img = (REAL*) mxGetData(Img); + if (numdim != 3) + { + mexPrintf("input has to be 3dim\n"); + return; + } + int w = dims[0]; + int h = dims[1]; + int d = dims[2]; + + const mxArray *BG; + BG = prhs[pcnt++]; + const int *dimsbg = mxGetDimensions(BG); + if (dimsbg[0] != 6) + { + mexPrintf("blurrguide has to be rank 2\n"); + return; + } + REAL *blurrg = (REAL*) mxGetData(BG); + + + const mxArray *Factors; + Factors = prhs[pcnt++]; + REAL *factors = (REAL*) mxGetData(Factors); + REAL alpha = factors[0]; + int itmax = (int) factors[1]; + + plhs[0] = mxCreateNumericArray(3,dims,mxGetClassID(Img),mxREAL); + REAL *result = (REAL*) mxGetData(plhs[0]); + memcpy(result,img,sizeof(REAL)*w*h*d); + + + + int ndims[4]; + ndims[0] = 3; ndims[1] = dims[0]; ndims[2] = dims[1]; ndims[3] = dims[2]; + const mxArray *Tmp1 = mxCreateNumericArray(4,ndims,mxGetClassID(Img),mxREAL); + REAL *tmp1 = (REAL*) mxGetData(Tmp1); + + + + int wh = w*h; + + REAL *D = blurrg; + for (int k = 0; k < itmax ; k++) + { + + for (int z = 1; z < d-1; z++) + for (int y = 1; y < h-1; y++) + for (int x = 1; x < w-1; x++) + { + int idx = z*wh + y*w + x; + REAL gx = result[idx-1] - result[idx]; + REAL gy = result[idx-w] - result[idx]; + REAL gz = result[idx-wh] - result[idx]; + + REAL *tmp = &(tmp1[3*idx]); + tmp[0] = gx*D[0] + gy*D[3] + gz*D[4]; + tmp[1] = gx*D[3] + gy*D[1] + gz*D[5]; + tmp[2] = gx*D[4] + gy*D[5] + gz*D[2]; + } + + for (int z = 1; z < d-1; z++) + for (int y = 1; y < h-1; y++) + for (int x = 1; x < w-1; x++) + { + int idx = (z*wh + y*w + x); + REAL *tmp = &(tmp1[3*idx]); + result[idx] = result[idx] - alpha*(tmp[3]-tmp[0] + tmp[1+3*w]-tmp[1] + tmp[2+3*wh]-tmp[2]); + } + /* + + + for (int z = 1; z < d-1; z++) + for (int y = 1; y < h-1; y++) + for (int x = 1; x < w-1; x++) + { + int idx = z*wh + y*w + x; + REAL gxx = result[idx+1] + result[idx-1] - 2*result[idx]; + REAL gyy = result[idx+w] + result[idx-w] - 2*result[idx]; + REAL gzz = result[idx+wh] + result[idx-wh] - 2*result[idx]; + REAL gxy = (result[idx+1+w] - result[idx-1+w] - (result[idx+1-w] - result[idx-1-w]))/4; + REAL gxz = (result[idx+1+wh] - result[idx-1+wh] - (result[idx+1-wh] - result[idx-1-wh]))/4; + REAL gyz = (result[idx+w+wh] - result[idx-w+wh] - (result[idx+w-wh] - result[idx-w-wh]))/4; + + REAL *D = &(blurrg[6*idx]); + //tmp[0] = gxx*D[0] + gxy*D[3] + gxz*D[4]; + //tmp[1] = gxy*D[3] + gyy*D[1] + gyz*D[5]; + //tmp[2] = gxz*D[4] + gyz*D[5] + gzz*D[2]; + + result[idx] = result[idx] - alpha*( gxx*D[0] + gxy*D[3] + gxz*D[4] + gxy*D[3] + gyy*D[1] + gyz*D[5] + gxz*D[4] + gyz*D[5] + gzz*D[2]); + } + */ + } + + + + + + +} + + + + + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/auxilary_classes.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/auxilary_classes.cpp new file mode 100644 index 0000000000..a159e35324 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/auxilary_classes.cpp @@ -0,0 +1,638 @@ + +#ifndef _AUXCLASS +#define _AUXCLASS + +#ifndef INFINITY +#define INFINITY 1000000000 +#endif + +//#define INLINE __attribute__ ((always_inline)) +#define INLINE inline + +//#define __SSE2 + + + +#ifdef __SSE2 +#include + +class pVector +{ + +private: + __m128 r; + +public: + + static float store[4]; + + pVector() + { + + } + + pVector(__m128 q) + { + r = q; + } + + + pVector(float x,float y,float z) + { + r = _mm_set_ps(0,z,y,x); + } + + + INLINE void storeXYZ() + { + _mm_store_ps(store,r); + } + + INLINE void setXYZ(float sx,float sy,float sz) + { + r = _mm_set_ps(0,sz,sy,sx); + } + + + + + INLINE void rand(float w,float h,float d) + { + float x = mtrand.frand()*w; + float y = mtrand.frand()*h; + float z = mtrand.frand()*d; + r = _mm_set_ps(0,z,y,x); + + } + + INLINE void rand_sphere() + { + r = _mm_set_ps(0,mtrand.frandn(),mtrand.frandn(),mtrand.frandn()); + normalize(); + } + + INLINE void normalize() + { + __m128 q = _mm_mul_ps(r,r); + _mm_store_ps(store,q); + float norm = sqrt(store[0]+store[1]+store[2]) + 0.00000001; + q = _mm_set_ps1(1/norm); + r = _mm_mul_ps(r,q); + + } + + INLINE float norm_square() + { + __m128 q = _mm_mul_ps(r,r); + _mm_store_ps(store,q); + return store[0]+store[1]+store[2]; + } + + INLINE void distortn(float sigma) + { + __m128 s = _mm_set_ps(0,mtrand.frandn(),mtrand.frandn(),mtrand.frandn()); + __m128 q = _mm_set_ps1(sigma); + r = _mm_add_ps(r,_mm_mul_ps(q,s)); + } + + INLINE pVector operator*(float s) + { + __m128 q = _mm_set_ps1(s); + return pVector(_mm_mul_ps(q,r)); + } + + INLINE void operator*=(float &s) + { + __m128 q = _mm_set_ps1(s); + r = _mm_mul_ps(q,r); + } + + INLINE pVector operator+(pVector R) + { + return pVector(_mm_add_ps(R.r,r)); + } + + INLINE void operator+=(pVector R) + { + r = _mm_add_ps(r,R.r); + } + + INLINE pVector operator-(pVector R) + { + return pVector(_mm_sub_ps(r,R.r)); + } + + INLINE void operator-=(pVector R) + { + r = _mm_sub_ps(r,R.r); + } + + INLINE pVector operator/(float &s) + { + __m128 q = _mm_set_ps1(s); + return pVector(_mm_div_ps(q,r)); + } + + INLINE void operator/=(float &s) + { + __m128 q = _mm_set_ps1(s); + r = _mm_div_ps(r,q); + } + + INLINE float operator*(pVector R) + { + __m128 q = _mm_mul_ps(r,R.r); + _mm_store_ps(store,q); + return store[0]+store[1]+store[2]; + } + + +}; + +float pVector::store[4]; + +#else + +class pVector +{ +private: + float x; + float y; + float z; + +public: + + pVector() + { + + } + + pVector(float x,float y,float z) + { + this->x = x; + this->y = y; + this->z = z; + } + + INLINE void SetXYZ(float sx,float sy, float sz) + { + x = sx; + y = sy; + z = sz; + } + + INLINE float* GetXYZ() + { + float xyz[3] = {x,y,z}; + return xyz; + } + + INLINE float GetX() + { + return x; + } + + INLINE float GetY() + { + return y; + } + + INLINE float GetZ() + { + return z; + } + + INLINE void rand(float w, float h, float d) + { + this->x = mtrand.frand()*w; + this->y = mtrand.frand()*h; + this->z = mtrand.frand()*d; + } + + INLINE void rand_sphere() + { + this->x = mtrand.frandn(); + this->y = mtrand.frandn(); + this->z = mtrand.frandn(); + normalize(); + } + + INLINE void normalize() + { + float norm = sqrt(x*x+y*y+z*z)+ 0.00000001; + *this /= norm; + } + + INLINE float norm_square() + { + return x*x + y*y + z*z; + } + + INLINE void distortn(float sigma) + { + x += sigma*mtrand.frandn(); + y += sigma*mtrand.frandn(); + z += sigma*mtrand.frandn(); + } + + INLINE float operator[](int index) + { + switch(index) + { + case 0: + return x; + case 1: + return y; + case 2: + return z; + } + } + + INLINE pVector operator*(float s) + { + return pVector(s*x,s*y,s*z); + } + + INLINE void operator*=(float &s) + { + x *= s; + y *= s; + z *= s; + } + + INLINE pVector operator+(pVector R) + { + return pVector(x+R.x,y+R.y,z+R.z); + } + + INLINE void operator+=(pVector R) + { + x += R.x; + y += R.y; + z += R.z; + } + + INLINE pVector operator-(pVector R) + { + return pVector(x-R.x,y-R.y,z-R.z); + } + + INLINE void operator-=(pVector R) + { + x -= R.x; + y -= R.y; + z -= R.z; + } + + INLINE pVector operator/(float &s) + { + return pVector(x/s,y/s,z/s); + } + + INLINE void operator/=(float &s) + { + x /= s; + y /= s; + z /= s; + } + + + INLINE float operator*(pVector R) + { + return x*R.x+y*R.y+z*R.z; + } +}; + +#endif + + +class Particle +{ +public: + + Particle() + { + label = 0; + pID = -1; + mID = -1; + } + + ~Particle() + { + } + + pVector R; + pVector N; + float cap; + float len; + + int gridindex; // index in the grid where it is living + int ID; + int pID; + int mID; + + int label; + int numerator; +}; + + +class EnergyGradient +{ +public: + pVector gR; + pVector gN; + + INLINE float norm2() + { + return gR.norm_square() + gN.norm_square(); + } +} ; + + +template +class SimpSamp +{ + + float *P; + int cnt; + +public: + T *objs; + + + SimpSamp() + { + P = (float*) malloc(sizeof(float)*1000); + objs = (T*) malloc(sizeof(T)*1000); + } + ~SimpSamp() + { + free(P); + free(objs); + } + + INLINE void clear() + { + cnt = 1; + P[0] = 0; + } + + INLINE void add(float p, T obj) + { + P[cnt] = P[cnt-1] + p; + objs[cnt-1] = obj; + cnt++; + } + + // INLINE int draw() + // { + // float r = mtrand.frand()*P[cnt-1]; + // for (int i = 1; i < cnt; i++) + // { + // if (r <= P[i]) + // return i-1; + // } + // return cnt-2; + // } + + INLINE int draw() + { + float r = mtrand.frand()*P[cnt-1]; + int j; + int rl = 1; + int rh = cnt-1; + while(rh != rl) + { + j = rl + (rh-rl)/2; + if (r < P[j]) + { + rh = j; + continue; + } + if (r > P[j]) + { + rl = j+1; + continue; + } + break; + } + return rh-1; + } + + + + + + INLINE T drawObj() + { + return objs[draw()]; + } + + INLINE bool isempty() + { + if (cnt == 1) + return true; + else + return false; + } + + + float probFor(int idx) + { + return (P[idx+1]-P[idx])/P[cnt-1]; + } + + float probFor(T &t) + { + for (int i = 1; i< cnt;i++) + { + if (t == objs[i-1]) + return probFor(i-1); + } + return 0; + } + + + +}; + + +class EndPoint +{ +public: + EndPoint() + {} + + EndPoint(Particle *p,int ep) + { + this->p = p; + this->ep = ep; + } + Particle *p; + int ep; + + inline bool operator==(EndPoint P) + { + return (P.p == p) && (P.ep == ep); + } +}; + +class Track +{ +public: + EndPoint track[1000]; + float energy; + float proposal_probability; + int length; + + void clear() + { + length = 0; + energy = 0; + proposal_probability = 1; + } + + + bool isequal(Track &t) + { + for (int i = 0; i < length;i++) + { + if (track[i].p != t.track[i].p || track[i].ep != t.track[i].ep) + return false; + } + return true; + } + +}; + +float getMax(float *arr, int cnt) +{ + float max = arr[0]; + for (int i = 1; i < cnt; i++) + { + if (arr[i] > max) + max = arr[i]; + } + return max; +} + + + +float getMin(float *arr, int cnt) +{ + float min = arr[0]; + for (int i = 1; i < cnt; i++) + { + if (arr[i] < min) + min = arr[i]; + } + return min; +} + + +int getArgMin(float *arr, int cnt) +{ + float min = arr[0]; + int idx = 0; + for (int i = 1; i < cnt; i++) + { + if (arr[i] < min) + { + min = arr[i]; + idx = i; + } + } + return idx; +} + + + + +inline float distLseg(pVector &R1,pVector &N1,pVector &R2,pVector &N2,float &len) +{ + + pVector D = R1-R2; + float beta = N1*N2; + float divisor = 1.001-beta*beta; + float gamma1 = N1*D; + float gamma2 = N2*D; + float t,u; + float EPdist[4]; + + pVector Q; + float dist = 102400000.0; + + while(true) + { + + t = -(gamma1+beta*gamma2) / divisor; + u = (gamma1*beta+gamma2) / divisor; + if (fabs(t) < len && fabs(u) < len) + { + Q = D + N1*t - N2*u; + dist = Q*Q; + break; + } + + beta = len*beta; + + t = beta - gamma1; + if (fabs(t) < len) + { + Q = D + N1*t - N2*len; + float d = Q*Q; + if (d < dist) dist = d; + } + + t = -beta - gamma1; + if (fabs(t) < len) + { + Q = D + N1*t + N2*len; + float d = Q*Q; + if (d < dist) dist = d; + } + + u = beta + gamma2; + if (fabs(u) < len) + { + Q = D + N1*len - N2*u; + float d = Q*Q; + if (d < dist) dist = d; + } + + u = -beta + gamma2; + if (fabs(u) < len) + { + Q = D - N1*len - N2*u; + float d = Q*Q; + if (d < dist) dist = d; + } + + if (dist != 102400000.0) + break; + + + EPdist[0] = beta + gamma1 - gamma2; + EPdist[1] = -beta + gamma1 + gamma2; + EPdist[2] = -beta - gamma1 - gamma2; + EPdist[3] = beta - gamma1 + gamma2; + int c = getArgMin(EPdist,4); + if (c==0) {t = +len; u = +len; } + if (c==1) {t = +len; u = -len; } + if (c==2) {t = -len; u = +len; } + if (c==3) {t = -len; u = -len; } + Q = D + N1*t - N2*u; + dist = Q*Q; + break; + + } + + + return dist; + +} + + +#endif + + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/pcRJMCMC.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/pcRJMCMC.cpp new file mode 100644 index 0000000000..f44e29260b --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/pcRJMCMC.cpp @@ -0,0 +1,279 @@ +#ifndef __pcRJMCMC_cpp__ +#define __pcRJMCMC_cpp__ + +//#include "mex.h" +//#include "matrix.h" +#include +#include +#include +#include + +using namespace std; + +#define float float +#define PI 3.1415926536 +//#define INFINITY 99999999999.0 + +//#define TIMING + +#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); +} + +Timing externalenergy_time; +Timing internalenergy_time; +Timing odfeval_time; +Timing total_time; + +Timing shiftproposal_time; +Timing birthproposal_time; +Timing deathproposal_time; +Timing capproposal_time; +Timing lenproposal_time; +Timing connproposal_time; + +PropStats deathstats; +PropStats birthstats; +PropStats connstats; +PropStats shiftstats; +PropStats capstats; +PropStats lenstats; + + +#endif + + + +#include "MersenneTwister.h" +MTRand mtrand; +float *BESSEL_APPROXCOEFF; + + +#include "EnergyComputer_connec.cpp" +//#include "EnergyComputer_center.cpp" +//#include "RJMCMC_singlegradprop.cpp" +#include "RJMCMC_randshift.cpp" + + +//void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) +//{ + +// if(nrhs != 7 && nrhs != 8) { + +// printf(" wrong usage!!!.\n\n"); +// return; +// } else if(nlhs>3) { +// printf("Too many output arguments\n"); +// return; +// } + +// int pcnt = 0; +// const mxArray *Points; +// Points = prhs[pcnt++]; +// int numPoints = mxGetN(Points); +// float *points = (float*) mxGetData(Points); + +// const mxArray *DataImg; +// DataImg = prhs[pcnt++]; +// float *dimg = (float*) mxGetData(DataImg); +// const int *dsize = mxGetDimensions(DataImg); + +// const mxArray *Mask; +// Mask = prhs[pcnt++]; +// float *mask = (float*) mxGetData(Mask); +// const int *dmsize = mxGetDimensions(Mask); +// int mask_oversamp_mult = dmsize[0]/dsize[1]; + + + +// const mxArray *VoxSize; +// VoxSize = prhs[pcnt++]; +// double *voxsize = (double*) mxGetPr(VoxSize); + + + +// const mxArray *Params; +// Params = prhs[pcnt++]; +// double *params = (double*) mxGetPr(Params); + +// float Temp = (float) params[0]; +// int numit = (int) params[1]; +// float conprob = (float) params[2]; +// float particle_weight = (float) params[3]; +// float particle_width = (float) params[4]; +// float particle_len = (float) params[5]; +// float chempot_connection = (float) params[6]; +// float chempot_particle = (float) params[7]; +// float inex_balance = (float) params[8]; +// float chempot2 = (float) params[9]; +// float meanval_sq = (float) params[10]; + +// const mxArray *BesselExpansion; +// BesselExpansion = prhs[pcnt++]; +// BESSEL_APPROXCOEFF = (float*) mxGetData(BesselExpansion); + + + +// // read spherical-interpolator data + +// const mxArray *sinterpstruct = prhs[pcnt++]; +// mxArray *Indices = mxGetField(sinterpstruct,0,"indices"); +// mxArray *BaryCoords = mxGetField(sinterpstruct,0,"barycoords"); +// mxArray *Beta = mxGetField(sinterpstruct,0,"beta"); +// mxArray *NumInterpPoints = mxGetField(sinterpstruct,0,"numpoints"); + +// float *indimg = (float*) mxGetData(Indices); +// const int *isize = mxGetDimensions(Indices); +// int totsz = isize[0]*isize[1]*isize[2]*isize[3]; +// int *indeximg = (int*) malloc(sizeof(int)*totsz); +// for (int k =0;k < totsz;k++) +// indeximg[k] = int(indimg[k])-1; +// float *barycoords = (float*) mxGetData(BaryCoords); +// float *beta = (float*) mxGetData(Beta); +// int nip = int(*((float*)mxGetData(NumInterpPoints))); + +// SphereInterpolator *sinterp = new SphereInterpolator(barycoords,indeximg,nip,isize[2],beta[0]); + +// double breakhandle = 0; +// const mxArray *BreakHandle; +// if (nrhs == 8) +// { +// BreakHandle = prhs[pcnt++]; +// breakhandle = *mxGetPr(BreakHandle); +// } + +// #ifdef TIMING +// externalenergy_time.clear(); +// internalenergy_time.clear(); +// odfeval_time.clear(); +// total_time.clear(); + +// shiftproposal_time.clear(); +// birthproposal_time.clear(); +// deathproposal_time.clear(); +// capproposal_time.clear(); +// lenproposal_time.clear(); +// connproposal_time.clear(); + +// deathstats.clear(); +// birthstats.clear(); +// connstats.clear(); +// shiftstats.clear(); +// capstats.clear(); +// lenstats.clear(); +// #endif + + + +// float cellsize = 2*particle_len; +// float curv_hardthres = 0.7; + +// fprintf(stderr,"setting up MH-sampler \n"); fflush(stderr); +// RJMCMC sampler(points,numPoints, dimg, dsize, voxsize, cellsize); +// fprintf(stderr,"setting up Energy-computer \n"); fflush(stderr); +// EnergyComputer encomp(dimg,dsize,voxsize,sinterp,&(sampler.pcontainer),mask,mask_oversamp_mult); + +// fprintf(stderr,"setting up parameters\n"); fflush(stderr); +// sampler.setParameters(Temp,numit,conprob,particle_len,curv_hardthres,chempot_particle); +// sampler.setEnergyComputer(&encomp); +// encomp.setParameters(particle_weight,particle_width,chempot_connection*particle_len*particle_len,particle_len,curv_hardthres,inex_balance,chempot2,meanval_sq); + +// fprintf(stderr,"starting to iterate\n"); fflush(stderr); +// sampler.iterate(breakhandle); + +// int cnt = sampler.pcontainer.pcnt; +// #ifdef TIMING +// mexPrintf("\nEnergy\n------------------------\n"); +// externalenergy_time.report("external "); +// odfeval_time.report("odfeval "); +// internalenergy_time.report("internal "); +// total_time.report_time("total energy comp. "); + +// mexPrintf("overhead for proposals und stuff:%.1fms\n", +// (total_time.time-(externalenergy_time.time+odfeval_time.time+internalenergy_time.time))/1000.0); + +// mexPrintf("\nProposals\n------------------------\n"); +// birthproposal_time.report("birth "); +// deathproposal_time.report("death "); +// shiftproposal_time.report("shift "); +// connproposal_time.report("conne "); +//// capproposal_time.report("capch "); +//// lenproposal_time.report("length "); +// mexPrintf("\n"); +// birthstats.report("birth "); +// deathstats.report("death "); +// shiftstats.report("shift "); +// connstats.report("conne "); +//// lenstats.report("length "); +//// capstats.report("capch "); +// mexPrintf("\n"); + + + +// #endif + + +// int dims[] = {sampler.attrcnt, sampler.pcontainer.pcnt}; +// plhs[0] = mxCreateNumericArray(2,dims,mxGetClassID(Points),mxfloat); +// float *npoints = (float*) mxGetData(plhs[0]); +// sampler.writeout(npoints); + + +// delete sinterp; +// free(indeximg); + + +//} + +#endif + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen.cpp new file mode 100644 index 0000000000..aaa49f676a --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen.cpp @@ -0,0 +1,124 @@ + +#include "mex.h" +#include "matrix.h" +#include +#include +#include +#include + +using namespace std; + +#define REAL float +#define PI 3.1415926536 + + + +#include "MersenneTwister.h" +//MTRand mtrand; +#include "auxilary_classes.cpp" + + +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; + } + + int ppcnt = 0; + const mxArray *Points; + Points = prhs[ppcnt++]; + int numPoints = mxGetN(Points); + REAL *points = (REAL*) mxGetData(Points); + + const mxArray *Param; + Param = prhs[ppcnt++]; + double *param = mxGetPr(Param); + + + Particle* particles = (Particle*) malloc(sizeof(Particle)*numPoints); + Particle particles_proposal[10000]; + + int i; + for (i = 0; i < numPoints; i++) + { + particles[i].R.setXYZ(points[3*i],points[3*i+1],points[3*i+2]); + } + pVector source = particles[0].R; + pVector sink = particles[numPoints-1].R; + + + int pcnt = numPoints; + REAL len = param[0]; // 0.9*minSpacing + REAL Leng = 0; + + particles_proposal[0].R = source; + REAL dtau = 0; + int cur_p = 1; + int cur_i = 1; + pVector dR; + REAL normdR; + for (;;) + { + while (dtau <= len && cur_p < pcnt) + { + dR = particles[cur_p].R - particles[cur_p-1].R; + normdR = sqrt(dR.norm_square()); + dtau += normdR; + Leng += normdR; + cur_p++; + } + + if (dtau >= len) + { + // proposal = current particle position - (current p - last p)*(current fibre length - len???)/(norm current p-last p) + particles_proposal[cur_i].R = particles[cur_p-1].R - dR*( (dtau-len)/normdR ); + } + else + { + particles_proposal[cur_i].R = sink; + break; + } + + dtau = dtau-len; + + cur_i++; + if (cur_i >= 10000) + { + mexPrintf("bugy"); + break; + } + + } + + + + + const int sz[] = {3, cur_i+1}; + + plhs[0] = mxCreateNumericArray(2,sz,mxGetClassID(Points),mxREAL); + REAL *pk = (REAL*) mxGetData(plhs[0]); + for (i = 0; i < cur_i+1 ; i++) + { + particles_proposal[i].R.storeXYZ(); + pk[3*i] = pVector::store[0]; + pk[3*i+1] = pVector::store[1]; + pk[3*i+2] = pVector::store[2]; + } + + + plhs[1] = mxCreateDoubleMatrix(1,1,mxREAL); + double *plen = mxGetPr(plhs[1]); + plen[0] = Leng; + + free(particles); +} + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen2.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen2.cpp new file mode 100644 index 0000000000..657f72491b --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen2.cpp @@ -0,0 +1,97 @@ +#include "mitkParticle.h" +#include "auxilary_classes.cpp" + +using namespace std; + +#define PI 3.1415926536 + +float squareNorm( vnl_vector_fixed v ) +{ + return v[0]*v[0]+v[1]*v[1]+v[2]*v[2]; +} + +mitk::Particle* createMitkParticle(Particle* p) +{ + mitk::Particle* particle = new mitk::Particle(); + pVector R = p->R; + vnl_vector_fixed pos; + pos[0] = R[0]-0.5; + pos[1] = R[1]-0.5; + pos[2] = R[2]-0.5; + particle->SetPosition(pos); + //MITK_INFO << "mitk: " << particle->GetPosition(); + pVector N = p->N; + vnl_vector_fixed dir; + dir[0] = N[0]; + dir[1] = N[1]; + dir[2] = N[2]; + particle->SetDirection(dir); + return particle; +} + +Particle* createParticle(Particle* particle){ + Particle* p = new Particle(); + p->R = pVector(particle->R); + p->N = pVector(particle->N); + return p; +} + +vector* ResampleFibers(vector* particleContainer, float len) +{ + //vector< mitk::Particle::Pointer >* outContainer = new::vector< mitk::Particle::Pointer >(); + vector< Particle* >* container = + new vector< Particle* >(); + + int numPoints = particleContainer->size(); + if(numPoints<2) + return container; + + Particle* source = createParticle(particleContainer->at(0)); + Particle* sink = createParticle(particleContainer->at(numPoints-1)); + + float Leng = 0; + + container->push_back(source); + // source->R.storeXYZ(); + // MITK_INFO << source->R.store[0] << ", " << source->R.store[1] << ", " << source->R.store[2]; + + float dtau = 0; + int cur_p = 1; + int cur_i = 1; + pVector dR; + float normdR; + + for (;;) + { + while (dtau <= len && cur_p < numPoints) + { + dR = particleContainer->at(cur_p)->R - particleContainer->at(cur_p-1)->R; + normdR = sqrt(dR.norm_square()); + dtau += normdR; + Leng += normdR; + cur_p++; + } + + if (dtau >= len) // if particles reach next voxel + { + // proposal = current particle position - (current p - last p)*(current fibre length - len???)/(norm current p-last p) + Particle* p = new Particle(); + p->R = particleContainer->at(cur_p-1)->R - dR*( (dtau-len)/normdR ); + p->N = dR; + // p->R.storeXYZ(); + // MITK_INFO << p->R.store[0] << ", " << p->R.store[1] << ", " << p->R.store[2]; + container->push_back(p); + } + else + { + container->push_back(sink); + break; + } + + dtau = dtau-len; + + cur_i++; + } + return container; +} + diff --git a/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen_new.cpp b/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen_new.cpp new file mode 100644 index 0000000000..3e98361a39 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GlobalTracking/reparametrize_arclen_new.cpp @@ -0,0 +1,69 @@ +#include + +using namespace std; + +#define PI 3.1415926536 + +float squareNorm( vnl_vector_fixed v ) +{ + return v[0]*v[0]+v[1]*v[1]+v[2]*v[2]; +} + +itk::VectorContainer::Pointer mexFunction(vector< mitk::Particle::Pointer >* particleContainer, float len) +{ + //vector< mitk::Particle::Pointer >* outContainer = new::vector< mitk::Particle::Pointer >(); + itk::VectorContainer::Pointer container = + itk::VectorContainer::New(); + + int numPoints = particleContainer->size(); + if(numPoints<2) + return container; + + mitk::Particle::Pointer source = particleContainer->at(0); + mitk::Particle::Pointer sink = particleContainer->at(numPoints-1); + + float Leng = 0; + + container->InsertElement(container->Size(), source); + + float dtau = 0; + int cur_p = 1; + int cur_i = 1; + vnl_vector_fixed dR; + float normdR; + + for (;;) + { + while (dtau <= len && cur_p < numPoints) + { + dR = particleContainer->at(cur_p)->GetPosition() - particleContainer->at(cur_p-1)->GetPosition(); // vector between two particle centers + normdR = squareNorm(dR); // square length of above vector + dtau += normdR; + Leng += normdR; + cur_p++; + } + + if (dtau >= len) // if particles reach next voxel + { + // proposal = current particle position - (current p - last p)*(current fibre length - len???)/(norm current p-last p) + //particles_proposal[cur_i].R = particles[cur_p-1].R - dR*( (dtau-len)/normdR ); + mitk::Particle::Pointer newSegment = mitk::Particle::New(); + newSegment->SetPosition(particleContainer->at(cur_p-1)->GetPosition() - dR*( (dtau-len)/normdR )); + dR.normalize(); + newSegment->SetDirection(dR); + newSegment->SetWidth(particleContainer->at(cur_p-1)->GetWidth()); + container->InsertElement(container->Size(), newSegment); + } + else + { + container->InsertElement(container->Size(), sink); + break; + } + + dtau = dtau-len; + + cur_i++; + } + return container; +} + diff --git a/Modules/DiffusionImaging/Tractography/itkGlobalTractographyFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGlobalTractographyFilter.cpp new file mode 100644 index 0000000000..ec4877cc62 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/itkGlobalTractographyFilter.cpp @@ -0,0 +1,472 @@ +#include "itkGlobalTractographyFilter.h" + +#include + +#include +#include "itkPointShell.h" + +#include "GlobalTracking/BuildFibres.cpp" + +#include + +#include +#include +#include +#include +#include +#include + +#include "GlobalTracking/reparametrize_arclen2.cpp" +#include + +struct LessDereference { + template + bool operator()(const T * lhs, const T * rhs) const { + return *lhs < *rhs; + } +}; + +namespace itk{ + + template< class TInputOdfImage, class TInputROIImage > + GlobalTractographyFilter< TInputOdfImage, TInputROIImage > + ::GlobalTractographyFilter(): + 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) + { + //this->m_MeasurementFrame.set_identity(); + this->SetNumberOfRequiredInputs(2); //Filter needs a DWI image + a Mask Image + } + + template< class TInputOdfImage, class TInputROIImage > + GlobalTractographyFilter< TInputOdfImage, TInputROIImage > + ::~GlobalTractographyFilter(){ + delete BESSEL_APPROXCOEFF; + if (m_Sampler!=NULL) + delete m_Sampler; + } + + template< class TInputOdfImage, class TInputROIImage > + void + GlobalTractographyFilter< 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 << "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 + GlobalTractographyFilter< TInputOdfImage, TInputROIImage > + ::BuildFibers(float* points, int numPoints) + { + MITK_INFO << "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); + + // initialize array of particles + CCAnalysis ccana(points, numPoints, spacing); + + // label the particles according to fiber affiliation and return number of fibers + int numFibers = ccana.iterate(m_FiberLength); + + if (numFibers<=0){ + MITK_INFO << "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 << "itkGlobalTractographyFilter: " << numFibers << " fibers accepted"; + } + + // fill output fiber bundle datastructure + template< class TInputOdfImage, class TInputROIImage > + typename GlobalTractographyFilter< TInputOdfImage, TInputROIImage >::FiberBundleType* + GlobalTractographyFilter< TInputOdfImage, TInputROIImage > + ::GetFiberBundle() + { + if (!m_AbortTracking) + { + m_BuildFibers = true; + while (m_BuildFibers){} + } + + return &m_FiberBundle; + } + + // get memory allocated for particle grid + template< class TInputOdfImage, class TInputROIImage > + float + GlobalTractographyFilter< TInputOdfImage, TInputROIImage > + ::GetMemoryUsage() + { + if (m_Sampler!=NULL) + return m_Sampler->m_ParticleGrid.GetMemoryUsage(); + return 0; + } + + // perform global tracking + template< class TInputOdfImage, class TInputROIImage > + void + GlobalTractographyFilter< 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 << "image size < 3 not supported"; + return; + } + + // 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 << "Image direction is not a rotation matrix. Tracking not possible!"; + return; + } + + // 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 + ifstream BaryCoords; + BaryCoords.open("FiberTrackingLUTBaryCoords.bin", 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 << "Unable to open barycoords file"; + return; + } + + ifstream Indices; + Indices.open("FiberTrackingLUTIndices.bin", 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 << "Unable to open indices file"; + return; + } + + // 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 << "not enough iterations!"; + return; + } + unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); + + // setup metropolis hastings sampler + MITK_INFO << "itkGlobalTractographyFilter: 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 << "itkGlobalTractographyFilter: 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; + + m_CurrentStep = step+1; + float temperature = m_TempStart * exp(alpha*(((1.0)*step)/((1.0)*m_Steps))); + MITK_INFO << "iterating step " << m_CurrentStep; + + m_Sampler->SetTemperature(temperature); + m_Sampler->Iterate(&m_ProposalAcceptance, &m_NumConnections, &m_NumParticles, &m_AbortTracking); + + MITK_INFO << "proposal acceptance: " << 100*m_ProposalAcceptance << "%"; + MITK_INFO << "particles: " << m_NumParticles; + MITK_INFO << "connections: " << m_NumConnections; + + 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 << "done generate data"; + } +} + + + + diff --git a/Modules/DiffusionImaging/Tractography/itkGlobalTractographyFilter.h b/Modules/DiffusionImaging/Tractography/itkGlobalTractographyFilter.h new file mode 100644 index 0000000000..45c31a9f09 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/itkGlobalTractographyFilter.h @@ -0,0 +1,164 @@ +#ifndef itkGlobalTractographyFilter_h +#define itkGlobalTractographyFilter_h + +#include "itkProcessObject.h" +#include "itkVectorContainer.h" +#include "itkImage.h" + +#include "GlobalTracking/pcRJMCMC.cpp" +#include "GlobalTracking/auxilary_classes.cpp" + +#include + +namespace itk{ + + template< class TInputQBallImage, class TInputROIImage > + class GlobalTractographyFilter : + public ProcessObject{ + public: + typedef GlobalTractographyFilter Self; + typedef ProcessObject Superclass; + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + itkNewMacro(Self); + itkTypeMacro( GlobalTractographyFilter, 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; + + itkSetMacro( TempStart, float ); + itkGetMacro( TempStart, float ); + + itkSetMacro( TempEnd, float ); + itkGetMacro( TempEnd, float ); + + itkSetMacro( NumIt, int ); + itkGetMacro( NumIt, int ); + + 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 ); + + itkGetMacro( Steps, 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); + + itkGetMacro(NumParticles, unsigned long); + itkGetMacro(NumConnections, unsigned long); + itkGetMacro(NumAcceptedFibers, int); + itkGetMacro(ProposalAcceptance, float); + + /** 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(); + float GetMemoryUsage(); + + protected: + + GlobalTractographyFilter(); + virtual ~GlobalTractographyFilter(); + + void ComputeFiberCorrelation(); + + void BuildFibers(float* points, int numPoints); + + // Input Images + typename InputQBallImageType::Pointer m_ItkQBallImage; + typename MaskImageType::Pointer m_MaskImage; + + // 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; + bool m_BuildFibers; + int m_Steps; + float m_Memory; + float m_ProposalAcceptance; + + RJMCMC* m_Sampler; + FiberBundleType m_FiberBundle; + unsigned long m_NumParticles; + unsigned long m_NumConnections; + }; +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkGlobalTractographyFilter.cpp" +#endif + +#endif diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 21a06e8a2b..59c190f5cc 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,95 +1,104 @@ -SET(CPP_FILES +SET(CPP_FILES - # Algorithms - Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h - Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h - Algorithms/itkTensorDerivedMeasurementsFilter.h - Algorithms/itkBrainMaskExtractionImageFilter.h - Algorithms/itkB0ImageExtractionImageFilter.h - Algorithms/itkTensorImageToDiffusionImageFilter.h - Algorithms/itkTensorToL2NormImageFilter.h - Algorithms/itkTractsToProbabilityImageFilter.h - # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundle IODataStructures/FiberBundle/mitkFiberBundle.cpp IODataStructures/FiberBundle/mitkFiberBundleWriter.cpp IODataStructures/FiberBundle/mitkFiberBundleReader.cpp IODataStructures/FiberBundle/mitkFiberBundleIOFactory.cpp IODataStructures/FiberBundle/mitkFiberBundleWriterFactory.cpp IODataStructures/FiberBundle/mitkFiberBundleSerializer.cpp + IODataStructures/FiberBundle/mitkParticle.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleMapper3D.cpp - + # Interactions Interactions/mitkFiberBundleInteractor.cpp ) -SET(H_FILES +SET(H_FILES + # Rendering Rendering/mitkDiffusionImageMapper.h + Rendering/mitkOdfVtkMapper2D.h + + # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h + + # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h - IODataStructures/FiberBundle/itkSlowPolyLineParametricPath.h - Rendering/mitkOdfVtkMapper2D.h + IODataStructures/FiberBundle/itkSlowPolyLineParametricPath.h + + # Tractography + Tractography/itkGlobalTractographyFilter.h + + # Algorithms + Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h + Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h + Algorithms/itkTensorDerivedMeasurementsFilter.h + Algorithms/itkBrainMaskExtractionImageFilter.h + Algorithms/itkB0ImageExtractionImageFilter.h + Algorithms/itkTensorImageToDiffusionImageFilter.h + Algorithms/itkTensorToL2NormImageFilter.h + Algorithms/itkTractsToProbabilityImageFilter.h ) SET( TOOL_FILES ) IF(WIN32) ENDIF(WIN32) -#MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp ) +#MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp )