diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui
index ce02a438eb..096e59cc5b 100644
--- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui
+++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui
@@ -1,216 +1,216 @@
QmitkPreprocessingViewControls
0
0
350
422
0
0
true
QmitkPreprocessingViewControls
-
Reduce Size
-
0
70
Multiple acquistions of one gradient direction can be averaged. Due to rounding errors, similar gradients often differ in the last decimal positions. The Merge radius allows to average them anyway by taking into account all directions within a certain radius.
true
-
QFrame::NoFrame
QFrame::Raised
0
-
Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured.
Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured.
Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured.
Merge radius
-
Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured.
Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured.
Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured.
6
2.000000000000000
0.000100000000000
0.001000000000000
-
false
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
Average redundant gradients
-
Non diffusion weighted image
-
0
30
Average and extract all images that were acquired without diffusion weighting.
true
-
false
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
Extract B0
-
Brain Mask
-
false
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
Estimate binary brain mask
-
Qt::Vertical
20
0
diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp
index b4e45203b0..fb4df74940 100644
--- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp
+++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionView.cpp
@@ -1,784 +1,779 @@
/*=========================================================================
Program: Medical Imaging & Interaction Toolkit
Module: $RCSfile$
Language: C++
Date: $Date: 2009-05-28 17:19:30 +0200 (Do, 28 Mai 2009) $
Version: $Revision: 17495 $
Copyright (c) German Cancer Research Center, Division of Medical and
Biological Informatics. All rights reserved.
See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
//#define MBILOG_ENABLE_DEBUG
#include "QmitkQBallReconstructionView.h"
#include "mitkDiffusionImagingConfigure.h"
// qt includes
#include
// itk includes
#include "itkTimeProbe.h"
// mitk includes
#include "mitkProgressBar.h"
#include "mitkStatusBar.h"
#include "mitkNodePredicateDataType.h"
#include "QmitkDataStorageComboBox.h"
#include "QmitkStdMultiWidget.h"
#include "itkDiffusionQballReconstructionImageFilter.h"
#include "itkAnalyticalDiffusionQballReconstructionImageFilter.h"
#include "itkVectorContainer.h"
#include "mitkQBallImage.h"
#include "mitkProperties.h"
#include "mitkVtkResliceInterpolationProperty.h"
#include "mitkLookupTable.h"
#include "mitkLookupTableProperty.h"
#include "mitkTransferFunction.h"
#include "mitkTransferFunctionProperty.h"
#include "mitkDataNodeObject.h"
#include "mitkOdfNormalizationMethodProperty.h"
#include "mitkOdfScaleByProperty.h"
#include "berryIStructuredSelection.h"
#include "berryIWorkbenchWindow.h"
#include "berryISelectionService.h"
#include
const std::string QmitkQBallReconstructionView::VIEW_ID =
"org.mitk.views.qballreconstruction";
#define DI_INFO MITK_INFO("DiffusionImaging")
typedef float TTensorPixelType;
const int QmitkQBallReconstructionView::nrconvkernels = 252;
using namespace berry;
struct QbrSelListener : ISelectionListener
{
berryObjectMacro(QbrSelListener);
QbrSelListener(QmitkQBallReconstructionView* view)
{
m_View = view;
}
void DoSelectionChanged(ISelection::ConstPointer selection)
{
// save current selection in member variable
m_View->m_CurrentSelection = selection.Cast();
// do something with the selected items
if(m_View->m_CurrentSelection)
{
bool foundDwiVolume = false;
// iterate selection
for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin();
i != m_View->m_CurrentSelection->End(); ++i)
{
// extract datatree node
if (mitk::DataNodeObject::Pointer nodeObj = i->Cast())
{
mitk::DataNode::Pointer node = nodeObj->GetDataNode();
// only look at interesting types
if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
{
foundDwiVolume = true;
}
}
}
m_View->m_Controls->m_ButtonStandard->setEnabled(foundDwiVolume);
}
}
void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection)
{
// check, if selection comes from datamanager
if (part)
{
QString partname(part->GetPartName().c_str());
if(partname.compare("Datamanager")==0)
{
// apply selection
DoSelectionChanged(selection);
}
}
}
QmitkQBallReconstructionView* m_View;
};
QmitkQBallReconstructionView::QmitkQBallReconstructionView()
: QmitkFunctionality(),
m_Controls(NULL),
m_MultiWidget(NULL)
{
}
QmitkQBallReconstructionView::QmitkQBallReconstructionView(const QmitkQBallReconstructionView& other)
{
Q_UNUSED(other);
throw std::runtime_error("Copy constructor not implemented");
}
//void QmitkQBallReconstructionView::OpactiyChanged(int value)
//{
// if (m_CurrentSelection)
// {
// if (mitk::DataNodeObject::Pointer nodeObj = m_CurrentSelection->Begin()->Cast())
// {
// mitk::DataNode::Pointer node = nodeObj->GetDataNode();
// if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
// {
// node->SetIntProperty("DisplayChannel", value);
// mitk::RenderingManager::GetInstance()->RequestUpdateAll();
// }
// }
// }
//}
//
//void QmitkQBallReconstructionView::OpactiyActionChanged()
//{
// if (m_CurrentSelection)
// {
// if (mitk::DataNodeObject::Pointer nodeObj = m_CurrentSelection->Begin()->Cast())
// {
// mitk::DataNode::Pointer node = nodeObj->GetDataNode();
// if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
// {
// int displayChannel = 0.0;
// if(node->GetIntProperty("DisplayChannel", displayChannel))
// {
// m_OpacitySlider->setValue(displayChannel);
// }
// }
// }
// }
//
// MITK_INFO << "changed";
//}
QmitkQBallReconstructionView::~QmitkQBallReconstructionView()
{
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->RemovePostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener);
}
void QmitkQBallReconstructionView::CreateQtPartControl(QWidget *parent)
{
if (!m_Controls)
{
// create GUI widgets
m_Controls = new Ui::QmitkQBallReconstructionViewControls;
m_Controls->setupUi(parent);
this->CreateConnections();
QStringList items;
items << "2" << "4" << "6" << "8";
m_Controls->m_QBallReconstructionMaxLLevelComboBox->addItems(items);
m_Controls->m_QBallReconstructionMaxLLevelComboBox->setCurrentIndex(1);
MethodChoosen(m_Controls->m_QBallReconstructionMethodComboBox->currentIndex());
- m_Controls->m_QBallReconstructionNumberThreadsSpinbox->setValue(8);
#ifndef DIFFUSION_IMAGING_EXTENDED
m_Controls->m_QBallReconstructionMethodComboBox->removeItem(3);
#endif
AdvancedCheckboxClicked();
// define data type for combobox
//m_Controls->m_ImageSelector->SetDataStorage( this->GetDefaultDataStorage() );
//m_Controls->m_ImageSelector->SetPredicate( mitk::NodePredicateDataType::New("DiffusionImage") );
}
m_SelListener = berry::ISelectionListener::Pointer(new QbrSelListener(this));
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->AddPostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener);
berry::ISelection::ConstPointer sel(
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager"));
m_CurrentSelection = sel.Cast();
m_SelListener.Cast()->DoSelectionChanged(sel);
}
void QmitkQBallReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget)
{
m_MultiWidget = &stdMultiWidget;
}
void QmitkQBallReconstructionView::StdMultiWidgetNotAvailable()
{
m_MultiWidget = NULL;
}
void QmitkQBallReconstructionView::CreateConnections()
{
if ( m_Controls )
{
connect( (QObject*)(m_Controls->m_ButtonStandard), SIGNAL(clicked()), this, SLOT(ReconstructStandard()) );
connect( (QObject*)(m_Controls->m_AdvancedCheckbox), SIGNAL(clicked()), this, SLOT(AdvancedCheckboxClicked()) );
connect( (QObject*)(m_Controls->m_QBallReconstructionMethodComboBox), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) );
}
}
void QmitkQBallReconstructionView::Activated()
{
QmitkFunctionality::Activated();
berry::ISelection::ConstPointer sel(
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager"));
m_CurrentSelection = sel.Cast();
m_SelListener.Cast()->DoSelectionChanged(sel);
}
void QmitkQBallReconstructionView::Deactivated()
{
QmitkFunctionality::Deactivated();
}
void QmitkQBallReconstructionView::ReconstructStandard()
{
int index = m_Controls->m_QBallReconstructionMethodComboBox->currentIndex();
#ifndef DIFFUSION_IMAGING_EXTENDED
if(index>=3)
{
index = index + 1;
}
#endif
switch(index)
{
case 0:
{
// Numerical
Reconstruct(0,0);
break;
}
case 1:
{
// Standard
Reconstruct(1,0);
break;
}
case 2:
{
// Solid Angle
Reconstruct(1,6);
break;
}
case 3:
{
// Constrained Solid Angle
Reconstruct(1,7);
break;
}
case 4:
{
// ADC
Reconstruct(1,4);
break;
}
case 5:
{
// Raw Signal
Reconstruct(1,5);
break;
}
}
}
void QmitkQBallReconstructionView::MethodChoosen(int method)
{
switch(method)
{
case 0:
m_Controls->m_Description->setText("Numerical recon. (Tuch2004)");
break;
case 1:
m_Controls->m_Description->setText("Spherical harmonics recon. (Descoteaux2007)");
break;
case 2:
m_Controls->m_Description->setText("SH recon. with solid angle consideration (Aganj2009)");
break;
case 3:
m_Controls->m_Description->setText("SH solid angle with non-neg. constraint (Goh2009)");
break;
case 4:
m_Controls->m_Description->setText("SH recon. of the plain ADC-profiles");
break;
case 5:
m_Controls->m_Description->setText("SH recon. of the raw diffusion signal");
break;
}
}
void QmitkQBallReconstructionView::AdvancedCheckboxClicked()
{
bool check = m_Controls->
m_AdvancedCheckbox->isChecked();
m_Controls->m_QBallReconstructionMaxLLevelTextLabel_2->setVisible(check);
m_Controls->m_QBallReconstructionMaxLLevelComboBox->setVisible(check);
m_Controls->m_QBallReconstructionLambdaTextLabel_2->setVisible(check);
m_Controls->m_QBallReconstructionLambdaLineEdit->setVisible(check);
m_Controls->m_QBallReconstructionThresholdLabel_2->setVisible(check);
m_Controls->m_QBallReconstructionThreasholdEdit->setVisible(check);
m_Controls->m_OutputB0Image->setVisible(check);
- m_Controls->m_QBallReconstructionNumberThreadsLabel_2->setVisible(check);
- m_Controls->m_QBallReconstructionNumberThreadsSpinbox->setVisible(check);
m_Controls->label_2->setVisible(check);
//m_Controls->textLabel1_2->setVisible(check);
//m_Controls->m_QBallReconstructionLambdaStepLineEdit->setVisible(check);
//m_Controls->textLabel1_3->setVisible(check);
m_Controls->frame_2->setVisible(check);
}
void QmitkQBallReconstructionView::Reconstruct(int method, int normalization)
{
if (m_CurrentSelection)
{
mitk::DataStorage::SetOfObjects::Pointer set =
mitk::DataStorage::SetOfObjects::New();
int at = 0;
for (IStructuredSelection::iterator i = m_CurrentSelection->Begin();
i != m_CurrentSelection->End();
++i)
{
if (mitk::DataNodeObject::Pointer nodeObj = i->Cast())
{
mitk::DataNode::Pointer node = nodeObj->GetDataNode();
if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
{
set->InsertElement(at++, node);
}
}
}
if(method == 0)
{
NumericalQBallReconstruction(set, normalization);
}
else
{
#if BOOST_VERSION / 100000 > 0
#if BOOST_VERSION / 100 % 1000 > 34
if(method == 1)
{
AnalyticalQBallReconstruction(set, normalization);
}
#else
std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
QMessageBox::warning(NULL,"ERROR","Boost 1.35 minimum required");
#endif
#else
std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
QMessageBox::warning(NULL,"ERROR","Boost 1.35 minimum required");
#endif
}
}
}
void QmitkQBallReconstructionView::NumericalQBallReconstruction
(mitk::DataStorage::SetOfObjects::Pointer inImages, int normalization)
{
try
{
itk::TimeProbe clock;
int nrFiles = inImages->size();
if (!nrFiles) return;
QString status;
mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles);
mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
std::vector nodes;
while ( itemiter != itemiterend ) // for all items
{
mitk::DiffusionImage* vols =
static_cast*>(
(*itemiter)->GetData());
std::string nodename;
(*itemiter)->GetStringProperty("name", nodename);
++itemiter;
// QBALL RECONSTRUCTION
clock.Start();
MBI_INFO << "QBall reconstruction ";
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf(
"QBall reconstruction for %s", nodename.c_str()).toAscii());
typedef itk::DiffusionQballReconstructionImageFilter
QballReconstructionImageFilterType;
QballReconstructionImageFilterType::Pointer filter =
QballReconstructionImageFilterType::New();
filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() );
- filter->SetNumberOfThreads( m_Controls->m_QBallReconstructionNumberThreadsSpinbox->value() );
filter->SetBValue(vols->GetB_Value());
filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->text().toFloat() );
switch(normalization)
{
case 0:
{
filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD);
break;
}
case 1:
{
filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO_B_VALUE);
break;
}
case 2:
{
filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_B_ZERO);
break;
}
case 3:
{
filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_NONE);
break;
}
default:
{
filter->SetNormalizationMethod(QballReconstructionImageFilterType::QBR_STANDARD);
}
}
filter->Update();
clock.Stop();
MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ;
// ODFs TO DATATREE
mitk::QBallImage::Pointer image = mitk::QBallImage::New();
image->InitializeByItk( filter->GetOutput() );
//image->SetImportVolume( filter->GetOutput()->GetBufferPointer(), 0, 0, mitk::Image::ImportMemoryManagementType::ManageMemory );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
mitk::DataNode::Pointer node=mitk::DataNode::New();
node->SetData( image );
QString newname;
newname = newname.append(nodename.c_str());
newname = newname.append("_QN%1").arg(normalization);
SetDefaultNodeProperties(node, newname.toStdString());
nodes.push_back(node);
// B-Zero TO DATATREE
if(m_Controls->m_OutputB0Image->isChecked())
{
mitk::Image::Pointer image4 = mitk::Image::New();
image4->InitializeByItk( filter->GetBZeroImage().GetPointer() );
image4->SetVolume( filter->GetBZeroImage()->GetBufferPointer() );
mitk::DataNode::Pointer node4=mitk::DataNode::New();
node4->SetData( image4 );
node4->SetProperty( "name", mitk::StringProperty::New(
QString(nodename.c_str()).append("_b0").toStdString()) );
nodes.push_back(node4);
}
mitk::ProgressBar::GetInstance()->Progress();
}
std::vector::iterator nodeIt;
for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt)
GetDefaultDataStorage()->Add(*nodeIt);
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii());
m_MultiWidget->RequestUpdate();
}
catch (itk::ExceptionObject &ex)
{
MBI_INFO << ex ;
return ;
}
}
void QmitkQBallReconstructionView::AnalyticalQBallReconstruction(
mitk::DataStorage::SetOfObjects::Pointer inImages,
int normalization)
{
try
{
itk::TimeProbe clock;
int nrFiles = inImages->size();
if (!nrFiles) return;
std::vector lambdas;
float minLambda = m_Controls->m_QBallReconstructionLambdaLineEdit->text().toFloat();
lambdas.push_back(minLambda);
int nLambdas = lambdas.size();
QString status;
mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles*nLambdas);
mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
std::vector* nodes
= new std::vector();
while ( itemiter != itemiterend ) // for all items
{
mitk::DiffusionImage* vols =
static_cast*>(
(*itemiter)->GetData());
std::string nodename;
(*itemiter)->GetStringProperty("name",nodename);
itemiter++;
// QBALL RECONSTRUCTION
clock.Start();
MBI_INFO << "QBall reconstruction ";
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf(
"QBall reconstruction for %s", nodename.c_str()).toAscii());
for(int i=0; im_QBallReconstructionMaxLLevelComboBox->currentIndex())
{
case 0:
{
TemplatedAnalyticalQBallReconstruction<2>(vols, currentLambda, nodename, nodes, normalization);
break;
}
case 1:
{
TemplatedAnalyticalQBallReconstruction<4>(vols, currentLambda, nodename, nodes, normalization);
break;
}
case 2:
{
TemplatedAnalyticalQBallReconstruction<6>(vols, currentLambda, nodename, nodes, normalization);
break;
}
case 3:
{
TemplatedAnalyticalQBallReconstruction<8>(vols, currentLambda, nodename, nodes, normalization);
break;
}
}
clock.Stop();
MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ;
mitk::ProgressBar::GetInstance()->Progress();
}
}
std::vector::iterator nodeIt;
for(nodeIt = nodes->begin(); nodeIt != nodes->end(); ++nodeIt)
GetDefaultDataStorage()->Add(*nodeIt);
m_MultiWidget->RequestUpdate();
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii());
}
catch (itk::ExceptionObject &ex)
{
MBI_INFO << ex ;
return ;
}
}
template
void QmitkQBallReconstructionView::TemplatedAnalyticalQBallReconstruction(
mitk::DiffusionImage* vols, float lambda,
std::string nodename, std::vector* nodes,
int normalization)
{
typedef itk::AnalyticalDiffusionQballReconstructionImageFilter
FilterType;
typename FilterType::Pointer filter = FilterType::New();
filter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() );
- filter->SetNumberOfThreads( m_Controls->m_QBallReconstructionNumberThreadsSpinbox->value() );
filter->SetBValue(vols->GetB_Value());
filter->SetThreshold( m_Controls->m_QBallReconstructionThreasholdEdit->text().toFloat() );
filter->SetLambda(lambda);
switch(normalization)
{
case 0:
{
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
break;
}
case 1:
{
filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO_B_VALUE);
break;
}
case 2:
{
filter->SetNormalizationMethod(FilterType::QBAR_B_ZERO);
break;
}
case 3:
{
filter->SetNormalizationMethod(FilterType::QBAR_NONE);
break;
}
case 4:
{
filter->SetNormalizationMethod(FilterType::QBAR_ADC_ONLY);
break;
}
case 5:
{
filter->SetNormalizationMethod(FilterType::QBAR_RAW_SIGNAL);
break;
}
case 6:
{
filter->SetNormalizationMethod(FilterType::QBAR_SOLID_ANGLE);
break;
}
case 7:
{
filter->SetNormalizationMethod(FilterType::QBAR_NONNEG_SOLID_ANGLE);
break;
}
default:
{
filter->SetNormalizationMethod(FilterType::QBAR_STANDARD);
}
}
filter->Update();
// ODFs TO DATATREE
mitk::QBallImage::Pointer image = mitk::QBallImage::New();
image->InitializeByItk( filter->GetOutput() );
image->SetVolume( filter->GetOutput()->GetBufferPointer() );
mitk::DataNode::Pointer node=mitk::DataNode::New();
node->SetData( image );
QString newname;
newname = newname.append(nodename.c_str());
newname = newname.append("_QA%1").arg(normalization);
SetDefaultNodeProperties(node, newname.toStdString());
nodes->push_back(node);
// mitk::Image::Pointer image5 = mitk::Image::New();
// image5->InitializeByItk( filter->GetODFSumImage().GetPointer() );
// image5->SetVolume( filter->GetODFSumImage()->GetBufferPointer() );
// mitk::DataNode::Pointer node5=mitk::DataNode::New();
// node5->SetData( image5 );
// node5->SetProperty( "name", mitk::StringProperty::New(
// QString(nodename.c_str()).append("_ODF").toStdString()) );
// nodes->push_back(node5);
// B-Zero TO DATATREE
if(m_Controls->m_OutputB0Image->isChecked())
{
mitk::Image::Pointer image4 = mitk::Image::New();
image4->InitializeByItk( filter->GetBZeroImage().GetPointer() );
image4->SetVolume( filter->GetBZeroImage()->GetBufferPointer() );
mitk::DataNode::Pointer node4=mitk::DataNode::New();
node4->SetData( image4 );
node4->SetProperty( "name", mitk::StringProperty::New(
QString(nodename.c_str()).append("_b0").toStdString()) );
nodes->push_back(node4);
}
}
void QmitkQBallReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name)
{
node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) );
node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) );
node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New());
node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New());
node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2));
node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1));
node->SetProperty( "visible", mitk::BoolProperty::New( true ) );
node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) );
node->SetProperty ("layer", mitk::IntProperty::New(100));
node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) );
//node->SetProperty( "opacity", mitk::FloatProperty::New(1.0f) );
node->SetProperty( "name", mitk::StringProperty::New(name) );
}
//node->SetProperty( "volumerendering", mitk::BoolProperty::New( false ) );
//node->SetProperty( "use color", mitk::BoolProperty::New( true ) );
//node->SetProperty( "texture interpolation", mitk::BoolProperty::New( true ) );
//node->SetProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() );
//node->SetProperty( "layer", mitk::IntProperty::New(0));
//node->SetProperty( "in plane resample extent by geometry", mitk::BoolProperty::New( false ) );
//node->SetOpacity(1.0f);
//node->SetColor(1.0,1.0,1.0);
//node->SetVisibility(true);
//node->SetProperty( "IsQBallVolume", mitk::BoolProperty::New( true ) );
//mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New();
//mitk::LevelWindow levelwindow;
//// levelwindow.SetAuto( image );
//levWinProp->SetLevelWindow( levelwindow );
//node->GetPropertyList()->SetProperty( "levelwindow", levWinProp );
//// add a default rainbow lookup table for color mapping
//if(!node->GetProperty("LookupTable"))
//{
// mitk::LookupTable::Pointer mitkLut = mitk::LookupTable::New();
// vtkLookupTable* vtkLut = mitkLut->GetVtkLookupTable();
// vtkLut->SetHueRange(0.6667, 0.0);
// vtkLut->SetTableRange(0.0, 20.0);
// vtkLut->Build();
// mitk::LookupTableProperty::Pointer mitkLutProp = mitk::LookupTableProperty::New();
// mitkLutProp->SetLookupTable(mitkLut);
// node->SetProperty( "LookupTable", mitkLutProp );
//}
//if(!node->GetProperty("binary"))
// node->SetProperty( "binary", mitk::BoolProperty::New( false ) );
//// add a default transfer function
//mitk::TransferFunction::Pointer tf = mitk::TransferFunction::New();
//node->SetProperty ( "TransferFunction", mitk::TransferFunctionProperty::New ( tf.GetPointer() ) );
//// set foldername as string property
//mitk::StringProperty::Pointer nameProp = mitk::StringProperty::New( name );
//node->SetProperty( "name", nameProp );
diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui
index 1b49026201..ce29e09017 100644
--- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui
+++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui
@@ -1,256 +1,230 @@
QmitkQBallReconstructionViewControls
0
0
350
385
0
0
true
QmitkQBallReconstructionViewControls
-
Reconstruction
-
Advanced Settings
-
QFrame::StyledPanel
QFrame::Raised
QFormLayout::AllNonFixedFieldsGrow
-
true
B0 Threshold
false
-
true
0
-
true
Output B0-Image
-
-
-
- true
-
-
- Number of Threads
-
-
- Qt::AlignRight|Qt::AlignTrailing|Qt::AlignVCenter
-
-
- false
-
-
-
- -
-
-
- true
-
-
- 4
-
-
-
- -
true
Spherical Harmonics:
- -
+
-
true
Maximum l-Level
false
- -
+
-
true
-1
- -
+
-
true
Regularization Parameter
Lambda
false
- -
+
-
true
0.006
-
2
-
Numerical
-
Standard (SH)
-
Solid Angle (SH)
-
Constraint Solid Angle (SH)
-
ADC-Profile only (SH)
-
Raw Signal only (SH)
-
TextLabel
-
false
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
Start Reconstruction
-
Qt::Vertical
20
0
diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp
index 2a8b7bdc61..026df14322 100644
--- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp
+++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionView.cpp
@@ -1,850 +1,847 @@
/*=========================================================================
Program: Medical Imaging & Interaction Toolkit
Module: $RCSfile$
Language: C++
Date: $Date: 2009-05-28 17:19:30 +0200 (Do, 28 Mai 2009) $
Version: $Revision: 17495 $
Copyright (c) German Cancer Research Center, Division of Medical and
Biological Informatics. All rights reserved.
See MITKCopyright.txt or http://www.mitk.org/copyright.html for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#include "QmitkTensorReconstructionView.h"
#include "mitkDiffusionImagingConfigure.h"
// qt includes
#include
// itk includes
#include "itkTimeProbe.h"
//#include "itkTensor.h"
// mitk includes
#include "mitkProgressBar.h"
#include "mitkStatusBar.h"
#include "mitkNodePredicateDataType.h"
#include "QmitkDataStorageComboBox.h"
#include "QmitkStdMultiWidget.h"
#include "mitkTeemDiffusionTensor3DReconstructionImageFilter.h"
#include "itkDiffusionTensor3DReconstructionImageFilter.h"
#include "itkTensorImageToDiffusionImageFilter.h"
#include "itkPointShell.h"
#include "itkVector.h"
#include "mitkProperties.h"
#include "mitkDataNodeObject.h"
#include "mitkOdfNormalizationMethodProperty.h"
#include "mitkOdfScaleByProperty.h"
#include "mitkDiffusionImageMapper.h"
#include "berryIStructuredSelection.h"
#include "berryIWorkbenchWindow.h"
#include "berryISelectionService.h"
#include
const std::string QmitkTensorReconstructionView::VIEW_ID =
"org.mitk.views.tensorreconstruction";
#define DI_INFO MITK_INFO("DiffusionImaging")
typedef float TTensorPixelType;
using namespace berry;
struct TrSelListener : ISelectionListener
{
berryObjectMacro(TrSelListener);
TrSelListener(QmitkTensorReconstructionView* view)
{
m_View = view;
}
void DoSelectionChanged(ISelection::ConstPointer selection)
{
// if(!m_View->IsVisible())
// return;
// save current selection in member variable
m_View->m_CurrentSelection = selection.Cast();
// do something with the selected items
if(m_View->m_CurrentSelection)
{
bool foundDwiVolume = false;
bool foundTensorVolume = false;
// iterate selection
for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin();
i != m_View->m_CurrentSelection->End(); ++i)
{
// extract datatree node
if (mitk::DataNodeObject::Pointer nodeObj = i->Cast())
{
mitk::DataNode::Pointer node = nodeObj->GetDataNode();
// only look at interesting types
if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
{
foundDwiVolume = true;
}
// only look at interesting types
if(QString("TensorImage").compare(node->GetData()->GetNameOfClass())==0)
{
foundTensorVolume = true;
}
}
}
m_View->m_Controls->m_ItkReconstruction->setEnabled(foundDwiVolume);
m_View->m_Controls->m_TeemReconstruction->setEnabled(foundDwiVolume);
m_View->m_Controls->m_TensorsToDWIButton->setEnabled(foundTensorVolume);
m_View->m_Controls->m_TensorsToQbiButton->setEnabled(foundTensorVolume);
}
}
void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection)
{
// check, if selection comes from datamanager
if (part)
{
QString partname(part->GetPartName().c_str());
if(partname.compare("Datamanager")==0)
{
// apply selection
DoSelectionChanged(selection);
}
}
}
QmitkTensorReconstructionView* m_View;
};
QmitkTensorReconstructionView::QmitkTensorReconstructionView()
: QmitkFunctionality(),
m_Controls(NULL),
m_MultiWidget(NULL)
{
}
QmitkTensorReconstructionView::QmitkTensorReconstructionView(const QmitkTensorReconstructionView& other)
{
Q_UNUSED(other)
throw std::runtime_error("Copy constructor not implemented");
}
QmitkTensorReconstructionView::~QmitkTensorReconstructionView()
{
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->RemovePostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener);
}
void QmitkTensorReconstructionView::CreateQtPartControl(QWidget *parent)
{
if (!m_Controls)
{
// create GUI widgets
m_Controls = new Ui::QmitkTensorReconstructionViewControls;
m_Controls->setupUi(parent);
this->CreateConnections();
- m_Controls->m_TensorReconstructionNumberThreadsSpinbox->setValue(8);
-
QStringList items;
items << "LLS (Linear Least Squares)"
<< "MLE (Maximum Likelihood)"
<< "NLS (Nonlinear Least Squares)"
<< "WLS (Weighted Least Squares)";
m_Controls->m_TensorEstimationTeemEstimationMethodCombo->addItems(items);
m_Controls->m_TensorEstimationTeemEstimationMethodCombo->setCurrentIndex(0);
m_Controls->m_TensorEstimationManualThreashold->setChecked(false);
m_Controls->m_TensorEstimationTeemSigmaEdit->setText("NaN");
m_Controls->m_TensorEstimationTeemNumItsSpin->setValue(1);
m_Controls->m_TensorEstimationTeemFuzzyEdit->setText("0.0");
m_Controls->m_TensorEstimationTeemMinValEdit->setText("1.0");
m_Controls->m_TensorEstimationTeemNumItsLabel_2->setEnabled(true);
m_Controls->m_TensorEstimationTeemNumItsSpin->setEnabled(true);
m_Controls->m_TensorsToDWIBValueEdit->setText("1000");
Advanced1CheckboxClicked();
Advanced2CheckboxClicked();
TeemCheckboxClicked();
#ifndef DIFFUSION_IMAGING_EXTENDED
m_Controls->m_TeemToggle->setVisible(false);
#endif
// define data type for combobox
//m_Controls->m_ImageSelector->SetDataStorage( this->GetDefaultDataStorage() );
//m_Controls->m_ImageSelector->SetPredicate( mitk::NodePredicateDataType::New("DiffusionImage") );
}
m_SelListener = berry::ISelectionListener::Pointer(new TrSelListener(this));
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->AddPostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener);
berry::ISelection::ConstPointer sel(
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager"));
m_CurrentSelection = sel.Cast();
m_SelListener.Cast()->DoSelectionChanged(sel);
}
void QmitkTensorReconstructionView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget)
{
berry::ISelection::ConstPointer sel(
this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager"));
m_CurrentSelection = sel.Cast();
m_SelListener.Cast()->DoSelectionChanged(sel);
m_MultiWidget = &stdMultiWidget;
}
void QmitkTensorReconstructionView::StdMultiWidgetNotAvailable()
{
m_MultiWidget = NULL;
}
void QmitkTensorReconstructionView::CreateConnections()
{
if ( m_Controls )
{
connect( (QObject*)(m_Controls->m_TeemToggle), SIGNAL(clicked()), this, SLOT(TeemCheckboxClicked()) );
connect( (QObject*)(m_Controls->m_ItkReconstruction), SIGNAL(clicked()), this, SLOT(ItkReconstruction()) );
connect( (QObject*)(m_Controls->m_TeemReconstruction), SIGNAL(clicked()), this, SLOT(TeemReconstruction()) );
connect( (QObject*)(m_Controls->m_TensorEstimationTeemEstimationMethodCombo), SIGNAL(currentIndexChanged(int)), this, SLOT(MethodChoosen(int)) );
connect( (QObject*)(m_Controls->m_Advanced1), SIGNAL(clicked()), this, SLOT(Advanced1CheckboxClicked()) );
connect( (QObject*)(m_Controls->m_Advanced2), SIGNAL(clicked()), this, SLOT(Advanced2CheckboxClicked()) );
connect( (QObject*)(m_Controls->m_TensorEstimationManualThreashold), SIGNAL(clicked()), this, SLOT(ManualThresholdClicked()) );
connect( (QObject*)(m_Controls->m_TensorsToDWIButton), SIGNAL(clicked()), this, SLOT(TensorsToDWI()) );
connect( (QObject*)(m_Controls->m_TensorsToQbiButton), SIGNAL(clicked()), this, SLOT(TensorsToQbi()) );
}
}
void QmitkTensorReconstructionView::TeemCheckboxClicked()
{
m_Controls->groupBox_3->setVisible(m_Controls->
m_TeemToggle->isChecked());
}
void QmitkTensorReconstructionView::Advanced1CheckboxClicked()
{
bool check = m_Controls->
m_Advanced1->isChecked();
m_Controls->frame->setVisible(check);
}
void QmitkTensorReconstructionView::Advanced2CheckboxClicked()
{
bool check = m_Controls->
m_Advanced2->isChecked();
m_Controls->frame_2->setVisible(check);
}
void QmitkTensorReconstructionView::ManualThresholdClicked()
{
m_Controls->m_TensorReconstructionThreasholdEdit_2->setEnabled(
m_Controls->m_TensorEstimationManualThreashold->isChecked());
}
void QmitkTensorReconstructionView::Activated()
{
QmitkFunctionality::Activated();
}
void QmitkTensorReconstructionView::Deactivated()
{
QmitkFunctionality::Deactivated();
}
void QmitkTensorReconstructionView::MethodChoosen(int method)
{
m_Controls->m_TensorEstimationTeemNumItsLabel_2->setEnabled(method==3);
m_Controls->m_TensorEstimationTeemNumItsSpin->setEnabled(method==3);
}
void QmitkTensorReconstructionView::ItkReconstruction()
{
Reconstruct(0);
}
void QmitkTensorReconstructionView::TeemReconstruction()
{
Reconstruct(1);
}
void QmitkTensorReconstructionView::Reconstruct(int method)
{
if (m_CurrentSelection)
{
mitk::DataStorage::SetOfObjects::Pointer set =
mitk::DataStorage::SetOfObjects::New();
int at = 0;
for (IStructuredSelection::iterator i = m_CurrentSelection->Begin();
i != m_CurrentSelection->End();
++i)
{
if (mitk::DataNodeObject::Pointer nodeObj = i->Cast())
{
mitk::DataNode::Pointer node = nodeObj->GetDataNode();
if(QString("DiffusionImage").compare(node->GetData()->GetNameOfClass())==0)
{
set->InsertElement(at++, node);
}
}
}
if(method == 0)
{
ItkTensorReconstruction(set);
}
if(method == 1)
{
TeemTensorReconstruction(set);
}
}
}
void QmitkTensorReconstructionView::ItkTensorReconstruction
(mitk::DataStorage::SetOfObjects::Pointer inImages)
{
try
{
itk::TimeProbe clock;
int nrFiles = inImages->size();
if (!nrFiles) return;
QString status;
mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles);
mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
std::vector nodes;
while ( itemiter != itemiterend ) // for all items
{
mitk::DiffusionImage* vols =
static_cast*>(
(*itemiter)->GetData());
std::string nodename;
(*itemiter)->GetStringProperty("name", nodename);
++itemiter;
// TENSOR RECONSTRUCTION
clock.Start();
MBI_INFO << "Tensor reconstruction ";
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf(
"Tensor reconstruction for %s", nodename.c_str()).toAscii());
typedef itk::DiffusionTensor3DReconstructionImageFilter<
DiffusionPixelType, DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType;
TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter =
TensorReconstructionImageFilterType::New();
tensorReconstructionFilter->SetGradientImage( vols->GetDirections(), vols->GetVectorImage() );
- tensorReconstructionFilter->SetNumberOfThreads( m_Controls->m_TensorReconstructionNumberThreadsSpinbox->value() );
tensorReconstructionFilter->SetBValue(vols->GetB_Value());
tensorReconstructionFilter->SetThreshold( m_Controls->m_TensorReconstructionThreasholdEdit->text().toFloat() );
tensorReconstructionFilter->Update();
clock.Stop();
MBI_DEBUG << "took " << clock.GetMeanTime() << "s.";
// TENSORS TO DATATREE
mitk::TensorImage::Pointer image = mitk::TensorImage::New();
typedef itk::Image, 3> TensorImageType;
TensorImageType::Pointer tensorImage;
tensorImage = tensorReconstructionFilter->GetOutput();
// Check the tensor for negative eigenvalues
if(m_Controls->m_CheckNegativeEigenvalues->isChecked())
{
typedef itk::ImageRegionIterator TensorImageIteratorType;
TensorImageIteratorType tensorIt(tensorImage, tensorImage->GetRequestedRegion());
tensorIt.GoToBegin();
while(!tensorIt.IsAtEnd())
{
typedef itk::DiffusionTensor3D TensorType;
//typedef itk::Tensor TensorType2;
TensorType tensor = tensorIt.Get();
// TensorType2 tensor2;
/*
for(int i=0; i SymEigenSystemType;
SymEigenSystemType eig (tensor2.GetVnlMatrix());
for(unsigned int i=0; iInitializeByItk( tensorImage.GetPointer() );
image->SetVolume( tensorReconstructionFilter->GetOutput()->GetBufferPointer() );
mitk::DataNode::Pointer node=mitk::DataNode::New();
node->SetData( image );
QString newname;
newname = newname.append(nodename.c_str());
newname = newname.append("_dti");
SetDefaultNodeProperties(node, newname.toStdString());
nodes.push_back(node);
mitk::ProgressBar::GetInstance()->Progress();
}
std::vector::iterator nodeIt;
for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt)
GetDefaultDataStorage()->Add(*nodeIt);
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii());
m_MultiWidget->RequestUpdate();
}
catch (itk::ExceptionObject &ex)
{
MBI_INFO << ex ;
return ;
}
}
void QmitkTensorReconstructionView::TeemTensorReconstruction
(mitk::DataStorage::SetOfObjects::Pointer inImages)
{
try
{
itk::TimeProbe clock;
int nrFiles = inImages->size();
if (!nrFiles) return;
QString status;
mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles);
mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
std::vector nodes;
while ( itemiter != itemiterend ) // for all items
{
mitk::DiffusionImage* vols =
static_cast*>(
(*itemiter)->GetData());
std::string nodename;
(*itemiter)->GetStringProperty("name", nodename);
++itemiter;
// TENSOR RECONSTRUCTION
clock.Start();
MBI_INFO << "Teem Tensor reconstruction ";
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf(
"Teem Tensor reconstruction for %s", nodename.c_str()).toAscii());
typedef mitk::TeemDiffusionTensor3DReconstructionImageFilter<
DiffusionPixelType, TTensorPixelType > TensorReconstructionImageFilterType;
TensorReconstructionImageFilterType::Pointer tensorReconstructionFilter =
TensorReconstructionImageFilterType::New();
tensorReconstructionFilter->SetInput( vols );
if(!m_Controls->m_TensorEstimationTeemSigmaEdit->text().contains(QString("NaN")))
tensorReconstructionFilter->SetSigma( m_Controls->m_TensorEstimationTeemSigmaEdit->text().toFloat() );
switch(m_Controls->m_TensorEstimationTeemEstimationMethodCombo->currentIndex())
{
// items << "LLS (Linear Least Squares)"
//<< "MLE (Maximum Likelihood)"
//<< "NLS (Nonlinear Least Squares)"
//<< "WLS (Weighted Least Squares)";
case 0:
tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsLLS);
break;
case 1:
tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsMLE);
break;
case 2:
tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsNLS);
break;
case 3:
tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsWLS);
break;
default:
tensorReconstructionFilter->SetEstimationMethod(mitk::TeemTensorEstimationMethodsLLS);
}
tensorReconstructionFilter->SetNumIterations( m_Controls->m_TensorEstimationTeemNumItsSpin->value() );
if(m_Controls->m_TensorEstimationManualThreashold->isChecked())
tensorReconstructionFilter->SetConfidenceThreshold( m_Controls->m_TensorReconstructionThreasholdEdit_2->text().toDouble() );
tensorReconstructionFilter->SetConfidenceFuzzyness( m_Controls->m_TensorEstimationTeemFuzzyEdit->text().toFloat() );
tensorReconstructionFilter->SetMinPlausibleValue( m_Controls->m_TensorEstimationTeemMinValEdit->text().toDouble() );
tensorReconstructionFilter->Update();
clock.Stop();
MBI_DEBUG << "took " << clock.GetMeanTime() << "s." ;
// TENSORS TO DATATREE
mitk::DataNode::Pointer node2=mitk::DataNode::New();
node2->SetData( tensorReconstructionFilter->GetOutputItk() );
QString newname;
newname = newname.append(nodename.c_str());
newname = newname.append("_dtix");
SetDefaultNodeProperties(node2, newname.toStdString());
nodes.push_back(node2);
mitk::ProgressBar::GetInstance()->Progress();
}
std::vector::iterator nodeIt;
for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt)
GetDefaultDataStorage()->Add(*nodeIt);
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii());
m_MultiWidget->RequestUpdate();
}
catch (itk::ExceptionObject &ex)
{
MBI_INFO << ex ;
return ;
}
}
void QmitkTensorReconstructionView::SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name)
{
node->SetProperty( "ShowMaxNumber", mitk::IntProperty::New( 500 ) );
node->SetProperty( "Scaling", mitk::FloatProperty::New( 1.0 ) );
node->SetProperty( "Normalization", mitk::OdfNormalizationMethodProperty::New());
node->SetProperty( "ScaleBy", mitk::OdfScaleByProperty::New());
node->SetProperty( "IndexParam1", mitk::FloatProperty::New(2));
node->SetProperty( "IndexParam2", mitk::FloatProperty::New(1));
node->SetProperty( "visible", mitk::BoolProperty::New( true ) );
node->SetProperty( "VisibleOdfs", mitk::BoolProperty::New( false ) );
node->SetProperty ("layer", mitk::IntProperty::New(100));
node->SetProperty( "DoRefresh", mitk::BoolProperty::New( true ) );
//node->SetProperty( "opacity", mitk::FloatProperty::New(1.0f) );
node->SetProperty( "name", mitk::StringProperty::New(name) );
}
//node->SetProperty( "volumerendering", mitk::BoolProperty::New( false ) );
//node->SetProperty( "use color", mitk::BoolProperty::New( true ) );
//node->SetProperty( "texture interpolation", mitk::BoolProperty::New( true ) );
//node->SetProperty( "reslice interpolation", mitk::VtkResliceInterpolationProperty::New() );
//node->SetProperty( "layer", mitk::IntProperty::New(0));
//node->SetProperty( "in plane resample extent by geometry", mitk::BoolProperty::New( false ) );
//node->SetOpacity(1.0f);
//node->SetColor(1.0,1.0,1.0);
//node->SetVisibility(true);
//node->SetProperty( "IsTensorVolume", mitk::BoolProperty::New( true ) );
//mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New();
//mitk::LevelWindow levelwindow;
//// levelwindow.SetAuto( image );
//levWinProp->SetLevelWindow( levelwindow );
//node->GetPropertyList()->SetProperty( "levelwindow", levWinProp );
//// add a default rainbow lookup table for color mapping
//if(!node->GetProperty("LookupTable"))
//{
// mitk::LookupTable::Pointer mitkLut = mitk::LookupTable::New();
// vtkLookupTable* vtkLut = mitkLut->GetVtkLookupTable();
// vtkLut->SetHueRange(0.6667, 0.0);
// vtkLut->SetTableRange(0.0, 20.0);
// vtkLut->Build();
// mitk::LookupTableProperty::Pointer mitkLutProp = mitk::LookupTableProperty::New();
// mitkLutProp->SetLookupTable(mitkLut);
// node->SetProperty( "LookupTable", mitkLutProp );
//}
//if(!node->GetProperty("binary"))
// node->SetProperty( "binary", mitk::BoolProperty::New( false ) );
//// add a default transfer function
//mitk::TransferFunction::Pointer tf = mitk::TransferFunction::New();
//node->SetProperty ( "TransferFunction", mitk::TransferFunctionProperty::New ( tf.GetPointer() ) );
//// set foldername as string property
//mitk::StringProperty::Pointer nameProp = mitk::StringProperty::New( name );
//node->SetProperty( "name", nameProp );
void QmitkTensorReconstructionView::TensorsToDWI()
{
if (m_CurrentSelection)
{
mitk::DataStorage::SetOfObjects::Pointer set =
mitk::DataStorage::SetOfObjects::New();
int at = 0;
for (IStructuredSelection::iterator i = m_CurrentSelection->Begin();
i != m_CurrentSelection->End();
++i)
{
if (mitk::DataNodeObject::Pointer nodeObj = i->Cast())
{
mitk::DataNode::Pointer node = nodeObj->GetDataNode();
if(QString("TensorImage").compare(node->GetData()->GetNameOfClass())==0)
{
set->InsertElement(at++, node);
}
}
}
DoTensorsToDWI(set);
}
}
void QmitkTensorReconstructionView::TensorsToQbi()
{
std::vector nodes = this->GetDataManagerSelection();
for (int i=0; i TensorPixelType;
typedef itk::Image< TensorPixelType, 3 > TensorImageType;
TensorImageType::Pointer itkvol = TensorImageType::New();
mitk::CastToItkImage(dynamic_cast(tensorImageNode->GetData()), itkvol);
typedef itk::TensorImageToQBallImageFilter< TTensorPixelType, TTensorPixelType > FilterType;
FilterType::Pointer filter = FilterType::New();
filter->SetInput( itkvol );
filter->Update();
typedef itk::Vector OutputPixelType;
typedef itk::Image OutputImageType;
mitk::QBallImage::Pointer image = mitk::QBallImage::New();
OutputImageType::Pointer outimg = filter->GetOutput();
image->InitializeByItk( outimg.GetPointer() );
image->SetVolume( outimg->GetBufferPointer() );
mitk::DataNode::Pointer node = mitk::DataNode::New();
node->SetData( image );
QString newname;
newname = newname.append(tensorImageNode->GetName().c_str());
newname = newname.append("_qbi");
node->SetName(newname.toAscii());
GetDefaultDataStorage()->Add(node);
}
}
void QmitkTensorReconstructionView::OnSelectionChanged( std::vector nodes )
{
if ( !this->IsVisible() )
return;
}
template
std::vector > QmitkTensorReconstructionView::MakeGradientList()
{
std::vector > retval;
vnl_matrix_fixed* U =
itk::PointShell >::DistributePointShell();
for(int i=0; i v;
v[0] = U->get(0,i); v[1] = U->get(1,i); v[2] = U->get(2,i);
retval.push_back(v);
}
return retval;
}
void QmitkTensorReconstructionView::DoTensorsToDWI
(mitk::DataStorage::SetOfObjects::Pointer inImages)
{
try
{
itk::TimeProbe clock;
int nrFiles = inImages->size();
if (!nrFiles) return;
QString status;
mitk::ProgressBar::GetInstance()->AddStepsToDo(nrFiles);
mitk::DataStorage::SetOfObjects::const_iterator itemiter( inImages->begin() );
mitk::DataStorage::SetOfObjects::const_iterator itemiterend( inImages->end() );
std::vector nodes;
while ( itemiter != itemiterend ) // for all items
{
std::string nodename;
(*itemiter)->GetStringProperty("name", nodename);
mitk::TensorImage* vol =
static_cast((*itemiter)->GetData());
++itemiter;
typedef float TTensorPixelType;
typedef itk::DiffusionTensor3D< TTensorPixelType > TensorPixelType;
typedef itk::Image< TensorPixelType, 3 > TensorImageType;
TensorImageType::Pointer itkvol = TensorImageType::New();
mitk::CastToItkImage(vol, itkvol);
typedef itk::TensorImageToDiffusionImageFilter<
TTensorPixelType, DiffusionPixelType > FilterType;
FilterType::GradientListType gradientList;
switch(m_Controls->m_TensorsToDWINumDirsSelect->currentIndex())
{
case 0:
gradientList = MakeGradientList<12>();
break;
case 1:
gradientList = MakeGradientList<42>();
break;
case 2:
gradientList = MakeGradientList<92>();
break;
case 3:
gradientList = MakeGradientList<162>();
break;
case 4:
gradientList = MakeGradientList<252>();
break;
case 5:
gradientList = MakeGradientList<362>();
break;
case 6:
gradientList = MakeGradientList<492>();
break;
case 7:
gradientList = MakeGradientList<642>();
break;
case 8:
gradientList = MakeGradientList<812>();
break;
case 9:
gradientList = MakeGradientList<1002>();
break;
default:
gradientList = MakeGradientList<92>();
}
double bVal = m_Controls->m_TensorsToDWIBValueEdit->text().toDouble();
// DWI ESTIMATION
clock.Start();
MBI_INFO << "DWI Estimation ";
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf(
"DWI Estimation for %s", nodename.c_str()).toAscii());
FilterType::Pointer filter = FilterType::New();
filter->SetInput( itkvol );
filter->SetBValue(bVal);
filter->SetGradientList(gradientList);
//filter->SetNumberOfThreads(1);
filter->Update();
clock.Stop();
MBI_DEBUG << "took " << clock.GetMeanTime() << "s.";
itk::Vector v;
v[0] = 0; v[1] = 0; v[2] = 0;
gradientList.push_back(v);
// TENSORS TO DATATREE
mitk::DiffusionImage::Pointer image = mitk::DiffusionImage::New();
image->SetVectorImage( filter->GetOutput() );
image->SetB_Value(bVal);
image->SetDirections(gradientList);
image->SetOriginalDirections(gradientList);
image->InitializeFromVectorImage();
mitk::DataNode::Pointer node=mitk::DataNode::New();
node->SetData( image );
mitk::DiffusionImageMapper::SetDefaultProperties(node);
QString newname;
newname = newname.append(nodename.c_str());
newname = newname.append("_dwi");
node->SetName(newname.toAscii());
nodes.push_back(node);
mitk::ProgressBar::GetInstance()->Progress();
}
std::vector::iterator nodeIt;
for(nodeIt = nodes.begin(); nodeIt != nodes.end(); ++nodeIt)
GetDefaultDataStorage()->Add(*nodeIt);
mitk::StatusBar::GetInstance()->DisplayText(status.sprintf("Finished Processing %d Files", nrFiles).toAscii());
m_MultiWidget->RequestUpdate();
}
catch (itk::ExceptionObject &ex)
{
MBI_INFO << ex ;
return ;
}
}
diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui
index 68d661b78e..7feacc3f22 100644
--- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui
+++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTensorReconstructionViewControls.ui
@@ -1,538 +1,524 @@
QmitkTensorReconstructionViewControls
0
0
345
836
0
0
true
QmitkTensorReconstructionViewControls
-
ITK Reconstruction
-
Advanced Settings
false
-
QFrame::StyledPanel
QFrame::Raised
+
+ QFormLayout::AllNonFixedFieldsGrow
+
-
-
-
- Number of Threads
-
-
- false
-
-
-
- -
-
-
- 4
-
-
-
- -
B0 Threshold
false
- -
+
-
0
- -
+
-
Check for negative eigenvalues
-
QFrame::StyledPanel
QFrame::Raised
-
false
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
ITK Tensor Reconstruction
-
Estimate Diffusion Image from Tensors
-
QFrame::StyledPanel
QFrame::Raised
QFormLayout::AllNonFixedFieldsGrow
6
6
9
-
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
B-Value
false
-
0
0
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
-
# Gradient Directions
-
3
-
12
-
42
-
92
-
162
-
252
-
362
-
492
-
642
-
812
-
1002
-
false
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
Diffusion Image Estimation
-
Estimate Q-Ball Image from Tensors
-
false
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+ Calculate ODF value as tensor value in the according direction
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
- sum of raw signal on equator, normalized to unit mass (Tuch et al. 2004)
+
Q-Ball Image Estimation
-
true
Teem Reconstruction
true
-
Teem Reconstruction
-
Advanced Settings
-
QFrame::StyledPanel
QFrame::Raised
QFormLayout::AllNonFixedFieldsGrow
-
important in case of method wls
# Iterations
false
-
-
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
Fuzzy confidence
false
-
0
0
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
how fuzzy the confidence boundary should be. By default, confidence boundary is perfectly sharp (float); default: "0"
-
minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0"
minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0"
minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0"
Min plausible value
false
-
Rician noise parameter (float)
Rician noise parameter (float)
Rician noise parameter (float)
Sigma
false
-
0
0
minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0"
minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0"
minimum plausible value (especially important for linear least squares estimation) (double); default: "1.0"
-
0
0
Rician noise parameter (float)
Rician noise parameter (float)
Rician noise parameter (float)
-
-
Method
-
B0-Threshold
-
false
0
-
false
Teem Tensor Reconstruction
-
Qt::Vertical
20
1150
diff --git a/Modules/DiffusionImaging/Algorithms/itkQBallToRgbImageFilter.h b/Modules/DiffusionImaging/Algorithms/itkQBallToRgbImageFilter.h
index 7f326d039b..5e093c135b 100644
--- a/Modules/DiffusionImaging/Algorithms/itkQBallToRgbImageFilter.h
+++ b/Modules/DiffusionImaging/Algorithms/itkQBallToRgbImageFilter.h
@@ -1,142 +1,145 @@
/*=========================================================================
Program: Medical Imaging & Interaction Toolkit
Language: C++
Date: $Date: 2009-07-14 19:11:20 +0200 (Tue, 14 Jul 2009) $
Version: $Revision: 18127 $
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 __itkQBallToRgbImageFilter_h
#define __itkQBallToRgbImageFilter_h
#include "itkUnaryFunctorImageFilter.h"
#include "itkOrientationDistributionFunction.h"
#include "itkRGBAPixel.h"
namespace itk
{
#define __IMG_DAT_ITEM__CEIL_ZERO_ONE__(val) (val) = \
( (val) < 0 ) ? ( 0 ) : ( ( (val)>1 ) ? ( 1 ) : ( (val) ) );
/** \class QBallToRgbImageFilter
*
*/
template ,3> >
class QBallToRgbImageFilter :
public ImageToImageFilter
{
public:
/** Standard class typedefs. */
typedef QBallToRgbImageFilter Self;
typedef ImageToImageFilter
Superclass;
typedef SmartPointer Pointer;
typedef SmartPointer ConstPointer;
typedef typename Superclass::InputImageType InputImageType;
typedef typename Superclass::OutputImageType OutputImageType;
typedef typename OutputImageType::PixelType OutputPixelType;
typedef typename TInputImage::PixelType InputPixelType;
typedef typename InputPixelType::ValueType InputValueType;
/** Run-time type information (and related methods). */
itkTypeMacro( QBallToRgbImageFilter, ImageToImageFilter );
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** Print internal ivars */
void PrintSelf(std::ostream& os, Indent indent) const
{ this->Superclass::PrintSelf( os, indent ); }
#ifdef ITK_USE_CONCEPT_CHECKING
/** Begin concept checking */
itkConceptMacro(InputHasNumericTraitsCheck,
(Concept::HasNumericTraits));
/** End concept checking */
#endif
protected:
QBallToRgbImageFilter() {};
virtual ~QBallToRgbImageFilter() {};
void GenerateData()
{
typename InputImageType::Pointer qballImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) );
typename OutputImageType::Pointer outputImage =
static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
typename InputImageType::RegionType region = qballImage->GetLargestPossibleRegion();
outputImage->SetSpacing( qballImage->GetSpacing() ); // Set the image spacing
outputImage->SetOrigin( qballImage->GetOrigin() ); // Set the image origin
outputImage->SetDirection( qballImage->GetDirection() ); // Set the image direction
outputImage->SetRegions( qballImage->GetLargestPossibleRegion());
outputImage->Allocate();
typedef ImageRegionConstIterator< InputImageType > QBallImageIteratorType;
QBallImageIteratorType qballIt(qballImage, qballImage->GetLargestPossibleRegion());
typedef ImageRegionIterator< OutputImageType > OutputImageIteratorType;
OutputImageIteratorType outputIt(outputImage, outputImage->GetLargestPossibleRegion());
qballIt.GoToBegin();
outputIt.GoToBegin();
while(!qballIt.IsAtEnd() && !outputIt.IsAtEnd()){
InputPixelType x = qballIt.Get();
typedef itk::OrientationDistributionFunction OdfType;
OdfType odf(x.GetDataPointer());
int pd = odf.GetPrincipleDiffusionDirection();
+ if (pd==-1)
+ MITK_ERROR << "ODF corrupted: GetPrincipleDiffusionDirection returned -1";
+
vnl_vector_fixed dir = OdfType::GetDirection(pd);
const float fa = odf.GetGeneralizedFractionalAnisotropy();
float r = fabs(dir[0]) * fa;
float g = fabs(dir[1]) * fa;
float b = fabs(dir[2]) * fa;
// float a = (fa-(m_OpacLevel-m_OpacWindow/2.0f))/m_OpacWindow;
float a = fa;
__IMG_DAT_ITEM__CEIL_ZERO_ONE__(r);
__IMG_DAT_ITEM__CEIL_ZERO_ONE__(g);
__IMG_DAT_ITEM__CEIL_ZERO_ONE__(b);
__IMG_DAT_ITEM__CEIL_ZERO_ONE__(a);
OutputPixelType out;
out.SetRed( r * 255.0f);
out.SetGreen( g * 255.0f);
out.SetBlue( b * 255.0f);
out.SetAlpha( a * 255.0f);
outputIt.Set(out);
++qballIt;
++outputIt;
}
}
private:
QBallToRgbImageFilter(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented
};
-
+
} // end namespace itk
-
+
#endif
diff --git a/Modules/DiffusionImaging/Algorithms/itkTensorImageToQBallImageFilter.txx b/Modules/DiffusionImaging/Algorithms/itkTensorImageToQBallImageFilter.txx
index bb2c8199dd..46697116ed 100644
--- a/Modules/DiffusionImaging/Algorithms/itkTensorImageToQBallImageFilter.txx
+++ b/Modules/DiffusionImaging/Algorithms/itkTensorImageToQBallImageFilter.txx
@@ -1,114 +1,115 @@
/*=========================================================================
Program: Tensor ToolKit - TTK
Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkTensorImageToQBallImageFilter.txx $
Language: C++
Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $
Version: $Revision: 68 $
Copyright (c) INRIA 2010. All rights reserved.
See LICENSE.txt for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#ifndef _itk_TensorImageToQBallImageFilter_txx_
#define _itk_TensorImageToQBallImageFilter_txx_
#endif
#include "itkTensorImageToQBallImageFilter.h"
#include
#include
#include
namespace itk
{
template
void
TensorImageToQBallImageFilter
::BeforeThreadedGenerateData()
{
typename OutputImageType::Pointer outImage = OutputImageType::New();
outImage->SetSpacing( this->GetInput()->GetSpacing() ); // Set the image spacing
outImage->SetOrigin( this->GetInput()->GetOrigin() ); // Set the image origin
outImage->SetDirection( this->GetInput()->GetDirection() ); // Set the image direction
outImage->SetLargestPossibleRegion( this->GetInput()->GetLargestPossibleRegion());
outImage->SetBufferedRegion( this->GetInput()->GetLargestPossibleRegion() );
outImage->SetRequestedRegion( this->GetInput()->GetLargestPossibleRegion() );
outImage->Allocate();
this->SetNumberOfRequiredOutputs (1);
this->SetNthOutput (0, outImage);
}
template
void
TensorImageToQBallImageFilter
::ThreadedGenerateData ( const OutputImageRegionType &outputRegionForThread, int threadId )
{
typedef ImageRegionIterator IteratorOutputType;
typedef ImageRegionConstIterator IteratorInputType;
unsigned long numPixels = outputRegionForThread.GetNumberOfPixels();
unsigned long step = numPixels/100;
unsigned long progress = 0;
IteratorOutputType itOut (this->GetOutput(0), outputRegionForThread);
IteratorInputType itIn (this->GetInput(), outputRegionForThread);
if( threadId==0 )
this->UpdateProgress (0.0);
while(!itIn.IsAtEnd())
{
if( this->GetAbortGenerateData() )
{
throw itk::ProcessAborted(__FILE__,__LINE__);
}
InputPixelType T = itIn.Get();
OutputPixelType out;
float tensorelems[6] = {
(float)T[0],
(float)T[1],
(float)T[2],
(float)T[3],
(float)T[4],
(float)T[5],
};
itk::DiffusionTensor3D tensor(tensorelems);
itk::OrientationDistributionFunction odf;
odf.InitFromTensor(tensor);
+ odf.Normalize();
for( unsigned int i=0; i0)
{
if( (progress%step)==0 )
{
this->UpdateProgress ( double(progress)/double(numPixels) );
}
}
++progress;
++itIn;
++itOut;
}
if( threadId==0 )
{
this->UpdateProgress (1.0);
}
MITK_INFO << "one thread finished Q-Ball estimation";
}
} // end of namespace
diff --git a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
index f668c26264..fda4d6c420 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
+++ b/Modules/DiffusionImaging/Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.cpp
@@ -1,899 +1,911 @@
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.txx,v $
Language: C++
Date: $Date: 2006-07-19 15:11:41 $
Version: $Revision: 1.11 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-This software is distributed WITHOUT ANY WARRANTY; without even
-the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+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 __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
#define __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
#include
#include
#include
#include
#include
#include
#include
#include
#include
#if BOOST_VERSION / 100000 > 0
#if BOOST_VERSION / 100 % 1000 > 34
#include
#endif
#endif
#include "itkPointShell.h"
namespace itk {
#define QBALL_ANAL_RECON_PI 3.14159265358979323846
template< class T, class TG, class TO, int L, int NODF>
AnalyticalDiffusionQballReconstructionImageFilter
::AnalyticalDiffusionQballReconstructionImageFilter() :
m_GradientDirectionContainer(NULL),
m_NumberOfGradientDirections(0),
m_NumberOfBaselineImages(1),
m_Threshold(NumericTraits< ReferencePixelType >::NonpositiveMin()),
m_BValue(1.0),
m_Lambda(0.0),
m_DirectionsDuplicated(false)
{
// At least 1 inputs is necessary for a vector image.
// For images added one at a time we need at least six
- this->SetNumberOfRequiredInputs( 1 );
+ this->SetNumberOfRequiredInputs( 1 );
}
- template<
- class TReferenceImagePixelType,
+ template<
+ class TReferenceImagePixelType,
class TGradientImagePixelType,
class TOdfPixelType,
int NOrderL,
int NrOdfDirections>
typename itk::AnalyticalDiffusionQballReconstructionImageFilter<
TReferenceImagePixelType,TGradientImagePixelType,TOdfPixelType,
- NOrderL,NrOdfDirections>::OdfPixelType
+ NOrderL,NrOdfDirections>::OdfPixelType
itk::AnalyticalDiffusionQballReconstructionImageFilter
-
- ::Normalize( OdfPixelType odf,
+ ::Normalize( OdfPixelType odf,
typename NumericTraits::AccumulateType b0 )
{
switch( m_NormalizationMethod )
{
case QBAR_STANDARD:
{
TOdfPixelType sum = 0;
for(int i=0; i0)
odf /= sum;
return odf;
break;
}
case QBAR_B_ZERO_B_VALUE:
{
for(int i=0; i0)
odf /= sum;
break;
}
case QBAR_NONNEG_SOLID_ANGLE:
{
break;
}
}
return odf;
}
- template<
- class TReferenceImagePixelType,
+ template<
+ class TReferenceImagePixelType,
class TGradientImagePixelType,
class TOdfPixelType,
int NOrderL,
int NrOdfDirections>
- vnl_vector
+ vnl_vector
itk::AnalyticalDiffusionQballReconstructionImageFilter
-
- ::PreNormalize( vnl_vector vec,
+ ::PreNormalize( vnl_vector vec,
typename NumericTraits::AccumulateType b0 )
{
switch( m_NormalizationMethod )
{
case QBAR_STANDARD:
{
return vec;
break;
}
case QBAR_B_ZERO_B_VALUE:
{
int n = vec.size();
for(int i=0; i= b0f)
- meas = b0f - 0.001;
- vec[i] = log(-log(meas/b0f));
+ if(vec[i] >= b0f)
+ vec[i] = b0f - 0.001;
+ vec[i] = log(-log(vec[i]/b0f));
}
return vec;
break;
}
}
return vec;
}
template< class T, class TG, class TO, int L, int NODF>
void AnalyticalDiffusionQballReconstructionImageFilter
::BeforeThreadedGenerateData()
{
- // If we have more than 2 inputs, then each input, except the first is a
+ // If we have more than 2 inputs, then each input, except the first is a
// gradient image. The number of gradient images must match the number of
// gradient directions.
//const unsigned int numberOfInputs = this->GetNumberOfInputs();
- // There need to be at least 6 gradient directions to be able to compute the
+ // There need to be at least 6 gradient directions to be able to compute the
// tensor basis
if( m_NumberOfGradientDirections < 6 )
{
itkExceptionMacro( << "At least 6 gradient directions are required" );
}
- // Input must be an itk::VectorImage.
+ // Input must be an itk::VectorImage.
std::string gradientImageClassName(
this->ProcessObject::GetInput(0)->GetNameOfClass());
if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 )
{
- itkExceptionMacro( <<
+ itkExceptionMacro( <<
"There is only one Gradient image. I expect that to be a VectorImage. "
<< "But its of type: " << gradientImageClassName );
}
this->ComputeReconstructionMatrix();
m_BZeroImage = BZeroImageType::New();
- typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
+ typename GradientImagesType::Pointer img = static_cast< GradientImagesType * >(
this->ProcessObject::GetInput(0) );
m_BZeroImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
m_BZeroImage->SetOrigin( img->GetOrigin() ); // Set the image origin
m_BZeroImage->SetDirection( img->GetDirection() ); // Set the image direction
m_BZeroImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
m_BZeroImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
m_BZeroImage->Allocate();
m_ODFSumImage = BZeroImageType::New();
m_ODFSumImage->SetSpacing( img->GetSpacing() ); // Set the image spacing
m_ODFSumImage->SetOrigin( img->GetOrigin() ); // Set the image origin
m_ODFSumImage->SetDirection( img->GetDirection() ); // Set the image direction
m_ODFSumImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion());
m_ODFSumImage->SetBufferedRegion( img->GetLargestPossibleRegion() );
m_ODFSumImage->Allocate();
if(m_NormalizationMethod == QBAR_SOLID_ANGLE ||
m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
{
m_Lambda = 0.0;
}
}
template< class T, class TG, class TO, int L, int NODF>
void AnalyticalDiffusionQballReconstructionImageFilter
::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread,
- int )
+ int )
{
- typename OutputImageType::Pointer outputImage =
+ typename OutputImageType::Pointer outputImage =
static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0));
ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread);
oit.GoToBegin();
ImageRegionIterator< BZeroImageType > oit2(m_BZeroImage, outputRegionForThread);
oit2.GoToBegin();
ImageRegionIterator< BlaImage > oit3(m_ODFSumImage, outputRegionForThread);
oit3.GoToBegin();
typedef ImageRegionConstIterator< GradientImagesType > GradientIteratorType;
typedef typename GradientImagesType::PixelType GradientVectorType;
typename GradientImagesType::Pointer gradientImagePointer = NULL;
// Would have liked a dynamic_cast here, but seems SGI doesn't like it
// The enum will ensure that an inappropriate cast is not done
- gradientImagePointer = static_cast< GradientImagesType * >(
+ gradientImagePointer = static_cast< GradientImagesType * >(
this->ProcessObject::GetInput(0) );
GradientIteratorType git(gradientImagePointer, outputRegionForThread );
git.GoToBegin();
// Compute the indicies of the baseline images and gradient images
std::vector baselineind; // contains the indicies of
// the baseline images
std::vector gradientind; // contains the indicies of
// the gradient images
for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
{
if(gdcit.Value().one_norm() <= 0.0)
{
baselineind.push_back(gdcit.Index());
}
else
{
gradientind.push_back(gdcit.Index());
}
}
if( m_DirectionsDuplicated )
{
int gradIndSize = gradientind.size();
for(int i=0; i::AccumulateType b0 = NumericTraits::Zero;
// Average the baseline image pixels
for(unsigned int i = 0; i < baselineind.size(); ++i)
{
b0 += b[baselineind[i]];
}
b0 /= this->m_NumberOfBaselineImages;
OdfPixelType odf(0.0);
vnl_vector B(m_NumberOfGradientDirections);
if( (b0 != 0) && (b0 >= m_Threshold) )
{
for( unsigned int i = 0; i< m_NumberOfGradientDirections; i++ )
{
B[i] = static_cast(b[gradientind[i]]);
}
B = PreNormalize(B, b0);
if(m_NormalizationMethod == QBAR_SOLID_ANGLE)
{
vnl_vector coeffs(m_NumberCoefficients);
coeffs = ( (*m_CoeffReconstructionMatrix) * B );
coeffs[0] += 1.0/(2.0*sqrt(QBALL_ANAL_RECON_PI));
odf = ( (*m_SphericalHarmonicBasisMatrix) * coeffs ).data_block();
}
else if(m_NormalizationMethod == QBAR_NONNEG_SOLID_ANGLE)
{
- /** this would be the place to implement a non-negative
+ /** this would be the place to implement a non-negative
* solver for quadratic programming problem:
* min .5*|| Bc-s ||^2 subject to -CLPc <= 4*pi*ones
* (refer to MICCAI 2009 Goh et al. "Estimating ODFs with PDF constraints")
* .5*|| Bc-s ||^2 == .5*c'B'Bc - x'B's + .5*s's
*/
itkExceptionMacro( << "Nonnegative Solid Angle not yet implemented");
// QuadProgPP::Matrix& G(m_G);
// int lenb = B.size();
// vnl_vector* s = new vnl_vector(lenb);
// for (int ii=0; ii g0_ = -1.0 * (*m_B_t) * (*s);
// QuadProgPP::Vector g0(g0_.data_block(),m_NumberCoefficients);
// try
// {
// QuadProgPP::QuadProg::solve_quadprog(G,g0,m_CE,m_ce0,m_CI,m_ci0,m_x);
// }
// catch(...)
// {
// m_x = 0;
// }
// vnl_vector coeffs(m_NumberCoefficients);
// for (int ii=0; ii
void AnalyticalDiffusionQballReconstructionImageFilter
::tofile2(vnl_matrix *pA, std::string fname)
{
vnl_matrix A = (*pA);
ofstream myfile;
std::locale C("C");
std::locale originalLocale = myfile.getloc();
myfile.imbue(C);
myfile.open (fname.c_str());
myfile << "A1=[";
for(int i=0; i
double AnalyticalDiffusionQballReconstructionImageFilter
::factorial(int number) {
if(number <= 1) return 1;
double result = 1.0;
for(int i=1; i<=number; i++)
result *= i;
return result;
}
template< class T, class TG, class TO, int L, int NODF>
void AnalyticalDiffusionQballReconstructionImageFilter
::Cart2Sph(double x, double y, double z, double *cart)
{
double phi, th, rad;
rad = sqrt(x*x+y*y+z*z);
th = atan2(z,sqrt(x*x+y*y));
phi = atan2(y,x);
th = -th+QBALL_ANAL_RECON_PI/2;
phi = -phi+QBALL_ANAL_RECON_PI;
cart[0] = phi;
cart[1] = th;
cart[2] = rad;
}
template< class T, class TG, class TO, int L, int NODF>
double AnalyticalDiffusionQballReconstructionImageFilter
::legendre0(int l)
{
if( l%2 != 0 )
{
return 0;
}
else
{
double prod1 = 1.0;
for(int i=1;i
double AnalyticalDiffusionQballReconstructionImageFilter
::spherical_harmonic(int m,int l,double theta,double phi, bool complexPart)
{
if( theta < 0 || theta > QBALL_ANAL_RECON_PI )
{
std::cout << "bad range" << std::endl;
return 0;
}
if( phi < 0 || phi > 2*QBALL_ANAL_RECON_PI )
{
std::cout << "bad range" << std::endl;
return 0;
}
double pml = 0;
double fac1 = factorial(l+m);
double fac2 = factorial(l-m);
if( m<0 )
{
#if BOOST_VERSION / 100000 > 0
#if BOOST_VERSION / 100 % 1000 > 34
pml = ::boost::math::legendre_p(l, -m, cos(theta));
#else
std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
#endif
#else
std::cout << "ERROR: Boost 1.35 minimum required" << std::endl;
#endif
double mypow = pow(-1.0,-m);
double myfac = (fac1/fac2);
pml *= mypow*myfac;
}
else
{
#if BOOST_VERSION / 100000 > 0
#if BOOST_VERSION / 100 % 1000 > 34
pml = ::boost::math::legendre_p(l, m, cos(theta));
#endif
#endif
}
//std::cout << "legendre(" << l << "," << m << "," << cos(theta) << ") = " << pml << std::endl;
double retval = sqrt(((2.0*(double)l+1.0)/(4.0*QBALL_ANAL_RECON_PI))*(fac2/fac1)) * pml;
if( !complexPart )
{
retval *= cos(m*phi);
}
else
{
retval *= sin(m*phi);
}
//std::cout << retval << std::endl;
return retval;
}
template< class T, class TG, class TO, int L, int NODF>
double AnalyticalDiffusionQballReconstructionImageFilter
::Yj(int m, int k, double theta, double phi)
{
if( -k <= m && m < 0 )
{
return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,false);
}
if( m == 0 )
return spherical_harmonic(0,k,theta,phi,false);
if( 0 < m && m <= k )
{
return sqrt(2.0) * spherical_harmonic(m,k,theta,phi,true);
}
return 0;
}
template< class T, class TG, class TO, int L, int NODF>
void AnalyticalDiffusionQballReconstructionImageFilter
::ComputeReconstructionMatrix()
{
//for(int i=-6;i<7;i++)
// std::cout << boost::math::legendre_p(6, i, 0.65657) << std::endl;
if( m_NumberOfGradientDirections < 6 )
{
itkExceptionMacro( << "Not enough gradient directions supplied. Need to supply at least 6" );
}
{
// check for duplicate diffusion gradients
bool warning = false;
for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
{
for(GradientDirectionContainerType::ConstIterator gdcit2 = this->m_GradientDirectionContainer->Begin();
gdcit2 != this->m_GradientDirectionContainer->End(); ++gdcit2)
{
if(gdcit1.Value() == gdcit2.Value() && gdcit1.Index() != gdcit2.Index())
{
itkWarningMacro( << "Some of the Diffusion Gradients equal each other. Corresponding image data should be averaged before calling this filter." );
warning = true;
break;
}
}
if (warning) break;
}
// handle acquisition schemes where only half of the spherical
// shell is sampled by the gradient directions. In this case,
// each gradient direction is duplicated in negative direction.
vnl_vector centerMass(3);
centerMass.fill(0.0);
int count = 0;
for(GradientDirectionContainerType::ConstIterator gdcit1 = this->m_GradientDirectionContainer->Begin();
gdcit1 != this->m_GradientDirectionContainer->End(); ++gdcit1)
{
if(gdcit1.Value().one_norm() > 0.0)
{
centerMass += gdcit1.Value();
count ++;
}
}
centerMass /= count;
if(centerMass.two_norm() > 0.1)
{
m_DirectionsDuplicated = true;
m_NumberOfGradientDirections *= 2;
}
}
vnl_matrix *Q
= new vnl_matrix(3, m_NumberOfGradientDirections);
{
int i = 0;
for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
{
if(gdcit.Value().one_norm() > 0.0)
{
double x = gdcit.Value().get(0);
double y = gdcit.Value().get(1);
double z = gdcit.Value().get(2);
double cart[3];
Cart2Sph(x,y,z,cart);
(*Q)(0,i) = cart[0];
(*Q)(1,i) = cart[1];
(*Q)(2,i++) = cart[2];
}
}
if(m_DirectionsDuplicated)
{
for(GradientDirectionContainerType::ConstIterator gdcit = this->m_GradientDirectionContainer->Begin();
gdcit != this->m_GradientDirectionContainer->End(); ++gdcit)
{
if(gdcit.Value().one_norm() > 0.0)
{
double x = gdcit.Value().get(0);
double y = gdcit.Value().get(1);
double z = gdcit.Value().get(2);
double cart[3];
Cart2Sph(x,y,z,cart);
(*Q)(0,i) = cart[0];
(*Q)(1,i) = cart[1];
(*Q)(2,i++) = cart[2];
}
}
}
}
int l = L;
m_NumberCoefficients = (int)(l*l + l + 2.0)/2.0 + l;
- vnl_matrix* B = new vnl_matrix(m_NumberOfGradientDirections,m_NumberCoefficients);
- vnl_matrix* _L = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
+ vnl_matrix* B = new vnl_matrix(m_NumberOfGradientDirections,m_NumberCoefficients);
+ vnl_matrix* _L = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
_L->fill(0.0);
- vnl_matrix* LL = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
+ vnl_matrix* LL = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
LL->fill(0.0);
- vnl_matrix* P = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
+ vnl_matrix* P = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
P->fill(0.0);
- vnl_matrix* Inv = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
+ vnl_matrix* Inv = new vnl_matrix(m_NumberCoefficients,m_NumberCoefficients);
P->fill(0.0);
vnl_vector* lj = new vnl_vector(m_NumberCoefficients);
m_LP = new vnl_vector(m_NumberCoefficients);
for(unsigned int i=0; i temp((*_L)*(*_L));
LL->update(*_L);
- *LL *= *_L;
+ *LL *= *_L;
//tofile2(LL,"LL");
for(int i=0; i(B->transpose());
//tofile2(&m_B_t,"m_B_t");
vnl_matrix B_t_B = (*m_B_t) * (*B);
//tofile2(&B_t_B,"B_t_B");
vnl_matrix lambdaLL(m_NumberCoefficients,m_NumberCoefficients);
lambdaLL.update((*LL));
lambdaLL *= m_Lambda;
//tofile2(&lambdaLL,"lLL");
vnl_matrix tmp( B_t_B + lambdaLL);
vnl_matrix_inverse *pseudoInverse
= new vnl_matrix_inverse( tmp );
(*Inv) = pseudoInverse->pinverse();
//tofile2(Inv,"Inv");
vnl_matrix temp((*Inv) * (*m_B_t));
double fac1 = (1.0/(16.0*QBALL_ANAL_RECON_PI*QBALL_ANAL_RECON_PI));
switch(m_NormalizationMethod)
{
case QBAR_ADC_ONLY:
case QBAR_RAW_SIGNAL:
break;
case QBAR_STANDARD:
case QBAR_B_ZERO_B_VALUE:
case QBAR_B_ZERO:
case QBAR_NONE:
temp = (*P) * temp;
break;
case QBAR_SOLID_ANGLE:
temp = fac1 * (*P) * (*_L) * temp;
break;
case QBAR_NONNEG_SOLID_ANGLE:
// m_G = QuadProgPP::Matrix(B_t_B.data_block(), B_t_B.rows(), B_t_B.cols());
// m_CE = QuadProgPP::Matrix((double)0,m_NumberCoefficients,0);
// m_ce0 = QuadProgPP::Vector((double)0,0);
// m_ci0 = QuadProgPP::Vector(4*QBALL_ANAL_RECON_PI, NODF);
// m_x = QuadProgPP::Vector(m_NumberCoefficients);
// (*m_LP) *= fac1;
break;
}
//tofile2(&temp,"A");
m_CoeffReconstructionMatrix = new vnl_matrix(m_NumberCoefficients,m_NumberOfGradientDirections);
for(int i=0; iodfs later
int NOdfDirections = NODF;
vnl_matrix_fixed* U =
itk::PointShell >::DistributePointShell();
m_SphericalHarmonicBasisMatrix = new vnl_matrix(NOdfDirections,m_NumberCoefficients);
- vnl_matrix* sphericalHarmonicBasisMatrix2
+ vnl_matrix* sphericalHarmonicBasisMatrix2
= new vnl_matrix(NOdfDirections,m_NumberCoefficients);
for(int i=0; i CI_t =
// (*sphericalHarmonicBasisMatrix2) * (*P) * (*_L);
// m_CI = QuadProgPP::Matrix(CI_t.transpose().data_block(), m_NumberCoefficients, NOdfDirections);
}
m_ReconstructionMatrix = new vnl_matrix(NOdfDirections,m_NumberOfGradientDirections);
*m_ReconstructionMatrix = (*m_SphericalHarmonicBasisMatrix) * (*m_CoeffReconstructionMatrix);
}
template< class T, class TG, class TO, int L, int NODF>
void AnalyticalDiffusionQballReconstructionImageFilter
- ::SetGradientImage( GradientDirectionContainerType *gradientDirection,
+ ::SetGradientImage( GradientDirectionContainerType *gradientDirection,
const GradientImagesType *gradientImage )
{
this->m_GradientDirectionContainer = gradientDirection;
unsigned int numImages = gradientDirection->Size();
this->m_NumberOfBaselineImages = 0;
for(GradientDirectionContainerType::Iterator it = this->m_GradientDirectionContainer->Begin();
it != this->m_GradientDirectionContainer->End(); it++)
{
if(it.Value().one_norm() <= 0.0)
{
this->m_NumberOfBaselineImages++;
}
else // Normalize non-zero gradient directions
{
it.Value() = it.Value() / it.Value().two_norm();
}
}
this->m_NumberOfGradientDirections = numImages - this->m_NumberOfBaselineImages;
- // ensure that the gradient image we received has as many components as
+ // ensure that the gradient image we received has as many components as
// the number of gradient directions
if( gradientImage->GetVectorLength() != this->m_NumberOfBaselineImages + m_NumberOfGradientDirections )
{
itkExceptionMacro( << m_NumberOfGradientDirections << " gradients + " << this->m_NumberOfBaselineImages
<< "baselines = " << m_NumberOfGradientDirections + this->m_NumberOfBaselineImages
<< " directions specified but image has " << gradientImage->GetVectorLength()
<< " components.");
}
- this->ProcessObject::SetNthInput( 0,
+ this->ProcessObject::SetNthInput( 0,
const_cast< GradientImagesType* >(gradientImage) );
}
template< class T, class TG, class TO, int L, int NODF>
void AnalyticalDiffusionQballReconstructionImageFilter
::PrintSelf(std::ostream& os, Indent indent) const
{
std::locale C("C");
std::locale originalLocale = os.getloc();
os.imbue(C);
Superclass::PrintSelf(os,indent);
os << indent << "OdfReconstructionMatrix: " << m_ReconstructionMatrix << std::endl;
if ( m_GradientDirectionContainer )
{
os << indent << "GradientDirectionContainer: "
<< m_GradientDirectionContainer << std::endl;
}
else
{
- os << indent <<
+ os << indent <<
"GradientDirectionContainer: (Gradient directions not set)" << std::endl;
}
- os << indent << "NumberOfGradientDirections: " <<
+ os << indent << "NumberOfGradientDirections: " <<
m_NumberOfGradientDirections << std::endl;
- os << indent << "NumberOfBaselineImages: " <<
+ os << indent << "NumberOfBaselineImages: " <<
m_NumberOfBaselineImages << std::endl;
os << indent << "Threshold for reference B0 image: " << m_Threshold << std::endl;
os << indent << "BValue: " << m_BValue << std::endl;
os.imbue( originalLocale );
}
}
#endif // __itkAnalyticalDiffusionQballReconstructionImageFilter_cpp
diff --git a/Modules/DiffusionImaging/Reconstruction/itkOrientationDistributionFunction.txx b/Modules/DiffusionImaging/Reconstruction/itkOrientationDistributionFunction.txx
index c084aee2d5..883e90b0be 100644
--- a/Modules/DiffusionImaging/Reconstruction/itkOrientationDistributionFunction.txx
+++ b/Modules/DiffusionImaging/Reconstruction/itkOrientationDistributionFunction.txx
@@ -1,1234 +1,1237 @@
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Language: C++
Date: $Date: 2008-03-10 22:48:13 $
Version: $Revision: 1.14 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
-This software is distributed WITHOUT ANY WARRANTY; without even
-the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+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 _itkOrientationDistributionFunction_txx
#define _itkOrientationDistributionFunction_txx
#include
#include
#include
#include
#include
#include
#include "itkPointShell.h"
//#include "itkNumericTraitsTensorPixel.h"
namespace itk
{
/*
-
+
#define INIT_STATIC_ODF_VARS(N_DIRS) \
_INIT_STATIC_ODF_VARS(float,N_DIRS) \
_INIT_STATIC_ODF_VARS(double,N_DIRS) \
-
+
#define _INIT_STATIC_ODF_VARS(PIXTYPE,N_DIRS) \
vtkPolyData* itk::OrientationDistributionFunction::m_BaseMesh = NULL; \
vnl_matrix_fixed* itk::OrientationDistributionFunction::m_Directions \
= itk::PointShell >::DistributePointShell(); \
std::vector< std::vector* >* itk::OrientationDistributionFunction::m_NeighborIdxs = NULL; \
std::vector* itk::OrientationDistributionFunction::m_HalfSphereIdxs = NULL; \
bool itk::OrientationDistributionFunction::m_Mutex = false; \
-
-
+
+
INIT_STATIC_ODF_VARS(40)
INIT_STATIC_ODF_VARS(60)
INIT_STATIC_ODF_VARS(80)
INIT_STATIC_ODF_VARS(100)
INIT_STATIC_ODF_VARS(150)
INIT_STATIC_ODF_VARS(200)
INIT_STATIC_ODF_VARS(250)
*/
template
vtkPolyData* itk::OrientationDistributionFunction::m_BaseMesh = NULL;
template
double itk::OrientationDistributionFunction::m_MaxChordLength = -1.0;
template
vnl_matrix_fixed* itk::OrientationDistributionFunction::m_Directions
= itk::PointShell >::DistributePointShell();
template
std::vector< std::vector* >* itk::OrientationDistributionFunction::m_NeighborIdxs = NULL;
template
std::vector< std::vector* >* itk::OrientationDistributionFunction::m_AngularRangeIdxs = NULL;
template
std::vector* itk::OrientationDistributionFunction::m_HalfSphereIdxs = NULL;
template
itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexBaseMesh;
template
itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexHalfSphereIdxs;
template
itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexNeighbors;
template
itk::SimpleFastMutexLock itk::OrientationDistributionFunction::m_MutexAngularRange;
#define ODF_PI 3.14159265358979323846
/*
* Assignment Operator
*/
template
OrientationDistributionFunction&
OrientationDistributionFunction
::operator= (const Self& r)
{
BaseArray::operator=(r);
return *this;
}
/*
* Assignment Operator from a scalar constant
*/
template
OrientationDistributionFunction&
OrientationDistributionFunction
::operator= (const ComponentType & r)
{
BaseArray::operator=(&r);
return *this;
}
/*
* Assigment from a plain array
*/
template
OrientationDistributionFunction&
OrientationDistributionFunction
::operator= (const ComponentArrayType r )
{
BaseArray::operator=(r);
return *this;
}
/**
* Returns a temporary copy of a vector
*/
template
- OrientationDistributionFunction
+ OrientationDistributionFunction
OrientationDistributionFunction
::operator+(const Self & r) const
{
Self result;
- for( unsigned int i=0; i
- OrientationDistributionFunction
+ OrientationDistributionFunction
OrientationDistributionFunction
::operator-(const Self & r) const
{
Self result;
- for( unsigned int i=0; i
- const OrientationDistributionFunction