diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake index 83c6bcfac2..725a26d759 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake @@ -1,92 +1,94 @@ SET(SRC_CPP_FILES QmitkODFDetailsWidget.cpp QmitkODFRenderWidget.cpp QmitkPartialVolumeAnalysisWidget.cpp + QmitkIVIMWidget.cpp ) SET(INTERNAL_CPP_FILES mitkPluginActivator.cpp QmitkQBallReconstructionView.cpp QmitkPreprocessingView.cpp QmitkDiffusionDicomImportView.cpp QmitkDiffusionQuantificationView.cpp QmitkTensorReconstructionView.cpp QmitkDiffusionImagingPublicPerspective.cpp QmitkControlVisualizationPropertiesView.cpp QmitkODFDetailsView.cpp QmitkGibbsTrackingView.cpp QmitkFiberBundleOperationsView.cpp QmitkFiberBundleDeveloperView.cpp QmitkPartialVolumeAnalysisView.cpp - + QmitkIVIMView.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/QmitkGibbsTrackingViewControls.ui src/internal/QmitkFiberBundleOperationsViewControls.ui src/internal/QmitkFiberBundleDeveloperViewControls.ui src/internal/QmitkPartialVolumeAnalysisViewControls.ui - + src/internal/QmitkIVIMViewControls.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/QmitkGibbsTrackingView.h src/internal/QmitkFiberBundleOperationsView.h src/internal/QmitkFiberBundleDeveloperView.h src/internal/QmitkPartialVolumeAnalysisView.h src/QmitkPartialVolumeAnalysisWidget.h - + src/internal/QmitkIVIMView.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/OdfDetails.png resources/GibbsTracking.png resources/FiberBundleOperations.png resources/PartialVolumeAnalysis_24.png + resources/IVIM_48.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 970a215860..45d0de0ebc 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml @@ -1,89 +1,96 @@ + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/IVIM_48.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/IVIM_48.png new file mode 100644 index 0000000000..de05f60c7c Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/IVIM_48.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc index f8e47d7a1f..36cf2a9503 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/QmitkDiffusionImaging.qrc @@ -1,31 +1,32 @@ qball.png tensor.png dwi.png dwiimport.png quantification.png reconodf.png recontensor.png texIntONIcon.png texIntOFFIcon.png vizControls.png Refresh_48.png QBallData24.png glyphsoff_C.png glyphsoff_S.png glyphsoff_T.png glyphson_C.png glyphson_S.png glyphson_T.png FiberBundle.png rectangle.png circle.png polygon.png color24.gif color48.gif color64.gif crosshair.png paint2.png + IVIM_48.png diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.cpp new file mode 100644 index 0000000000..4dc1f16b44 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.cpp @@ -0,0 +1,157 @@ +/*========================================================================= +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 "QmitkIVIMWidget.h" + +#include "mitkHistogramGenerator.h" + +#include +#include +#include +//#include + + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +//#include + +QmitkIVIMWidget::QmitkIVIMWidget( QWidget * parent ) + : QmitkPlotWidget(parent) +{ + // this->SetAxisTitle( QwtPlot::xBottom, "Grayvalue" ); + // this->SetAxisTitle( QwtPlot::yLeft, "Probability" ); + // this->Replot(); + m_Plot->setCanvasLineWidth(0); + m_Plot->setMargin(0); + + QwtLog10ScaleEngine* logScale = new QwtLog10ScaleEngine(); + m_Plot->setAxisScaleEngine(0, logScale); + + m_Plot->setAxisScale( 0, 0.15, 1.0 ); +} + + + +QmitkIVIMWidget::~QmitkIVIMWidget() +{ +} + + +void QmitkIVIMWidget::DrawGauss() +{ + +} + + +void QmitkIVIMWidget::ClearItemModel() +{ + +} + +std::vector QmitkIVIMWidget::vec(vnl_vector vector) +{ + std::vector retval(vector.size()); + for(unsigned int i=0; iClear(); + + QString s("f=%1, D=%2, D*=%3"); + s = s.arg(snap.currentF,4); + s = s.arg(snap.currentD,4); + s = s.arg(snap.currentDStar,4); + int curveId = this->InsertCurve( s.toAscii() ); + this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); + + curveId = this->InsertCurve( "ignored measurement points" ); + this->SetCurveData( curveId, vec(snap.bvalues), vec(snap.allmeas) ); + this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); + this->SetCurveSymbol(curveId, &QwtSymbol(QwtSymbol::Diamond, QColor(Qt::white), QColor(Qt::black), QSize(10,10))); + + if(snap.currentDStar != 0) + { + curveId = this->InsertCurve( "additional points second fit" ); + this->SetCurveData( curveId, vec(snap.bvals2), vec(snap.meas2) ); + this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); + this->SetCurveSymbol(curveId, &QwtSymbol(QwtSymbol::Diamond, QColor(Qt::black), QColor(Qt::black), QSize(10,10))); + } + + curveId = this->InsertCurve( "points first fit" ); + this->SetCurveData( curveId, vec(snap.bvals1), vec(snap.meas1) ); + this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); + this->SetCurveSymbol(curveId, &QwtSymbol(QwtSymbol::Diamond, QColor(Qt::red), QColor(Qt::red), QSize(10,10))); + + QPen pen; + pen.setColor( QColor(Qt::red) ); + pen.setWidth(2); + double maxb = snap.bvalues.max_value(); + vnl_vector xvals(2); + vnl_vector yvals(2); + xvals[0] = 0; + xvals[1] = maxb; + yvals[0] = 1-snap.currentFunceiled; + yvals[1] = yvals[0]*exp(-maxb * snap.currentD); + curveId = this->InsertCurve( "contribution of D to the signal" ); + this->SetCurveData( curveId, vec(xvals), vec(yvals) ); + this->SetCurvePen( curveId, pen ); + + if(snap.currentDStar != 0) + { + pen.setColor(Qt::black); + int nsampling = 50; + xvals.set_size(nsampling); + yvals.set_size(nsampling); + double f = 1-snap.currentFunceiled; + for(int i=0; iInsertCurve( "resulting fit of the model" ); + this->SetCurveData( curveId, vec(xvals), vec(yvals) ); + this->SetCurvePen( curveId, pen ); + } + +// QMargins margins; +// margins.setBottom(0); +// margins.setLeft(0); +// margins.setRight(0); +// margins.setTop(0); + + QwtLegend* legend = new QwtLegend(); +// legend->setContentsMargins(margins); + m_Plot->insertLegend(legend, QwtPlot::BottomLegend); + + this->Replot(); + +} + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.h new file mode 100644 index 0000000000..aa8ad3c104 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkIVIMWidget.h @@ -0,0 +1,73 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-05-15 18:09:46 +0200 (Fr, 15 Mai 2009) $ +Version: $Revision: 1.12 $ + +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 QmitkIVIMWidget_H_ +#define QmitkIVIMWidget_H_ + +#include "QmitkPlotWidget.h" + +#include "QmitkHistogram.h" +#include "QmitkExtExports.h" +#include "mitkImage.h" +#include "mitkPlanarFigure.h" + +#include +#include +#include + +#include + +#include +#include + +#include + +#include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" + + +/** + * \brief Widget for displaying image histograms based on the vtkQtChart + * framework + */ +class QmitkIVIMWidget : public QmitkPlotWidget +{ + +public: + + typedef itk::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter IVIMFilterType; + + typedef mitk::Image::HistogramType HistogramType; + typedef mitk::Image::HistogramType::ConstIterator HistogramConstIteratorType; + + void SetParameters( IVIMFilterType::IVIMSnapshot snap ); + + QmitkIVIMWidget( QWidget * /*parent = 0 */); + virtual ~QmitkIVIMWidget(); + + void DrawGauss(); + + void ClearItemModel(); + + std::vector< std::vector* > m_Vals; + +private: + + std::vector vec(vnl_vector vector); + +}; + +#endif /* QmitkIVIMWidget_H_ */ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.cpp new file mode 100644 index 0000000000..451c60bd3c --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.cpp @@ -0,0 +1,781 @@ +/*========================================================================= +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 "QmitkIVIMView.h" +#include "QmitkStdMultiWidget.h" + +// qt +#include "qmessagebox.h" +#include "qclipboard.h" + +// mitk +#include "mitkDiffusionImage.h" +#include "mitkImageCast.h" + +// itk +#include "itkScalarImageToHistogramGenerator.h" +#include "itkRegionOfInterestImageFilter.h" +#include "itkImageRegionConstIteratorWithIndex.h" + +// itk/mitk +#include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" +#include "itkRegularizedIVIMReconstructionFilter.h" +#include "mitkImageCast.h" + +const std::string QmitkIVIMView::VIEW_ID = "org.mitk.views.ivim"; + +QmitkIVIMView::QmitkIVIMView() + : QmitkFunctionality() + , m_Controls( 0 ) + , m_MultiWidget( NULL ) +{ +} + +QmitkIVIMView::~QmitkIVIMView() +{ +} + +void QmitkIVIMView::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::QmitkIVIMViewControls; + m_Controls->setupUi( parent ); + + connect( m_Controls->m_ButtonStart, SIGNAL(clicked()), this, SLOT(FittIVIMStart()) ); + connect( m_Controls->m_ButtonAutoThres, SIGNAL(clicked()), this, SLOT(AutoThreshold()) ); + + connect( m_Controls->m_MethodCombo, SIGNAL(currentIndexChanged(int)), this, SLOT(MethodCombo(int)) ); + + connect( m_Controls->m_DStarSlider, SIGNAL(valueChanged(int)), this, SLOT(DStarSlider(int)) ); + connect( m_Controls->m_BThreshSlider, SIGNAL(valueChanged(int)), this, SLOT(BThreshSlider(int)) ); + connect( m_Controls->m_S0ThreshSlider, SIGNAL(valueChanged(int)), this, SLOT(S0ThreshSlider(int)) ); + connect( m_Controls->m_NumItSlider, SIGNAL(valueChanged(int)), this, SLOT(NumItsSlider(int)) ); + connect( m_Controls->m_LambdaSlider, SIGNAL(valueChanged(int)), this, SLOT(LambdaSlider(int)) ); + + connect( m_Controls->m_DisplayResultsCheckbox, SIGNAL(clicked()), this, SLOT(Checkbox()) ); + connect( m_Controls->m_CheckDStar, SIGNAL(clicked()), this, SLOT(Checkbox()) ); + connect( m_Controls->m_CheckD, SIGNAL(clicked()), this, SLOT(Checkbox()) ); + connect( m_Controls->m_Checkf, SIGNAL(clicked()), this, SLOT(Checkbox()) ); + + connect( m_Controls->m_ChooseMethod, SIGNAL(clicked()), this, SLOT(ChooseMethod()) ); + connect( m_Controls->m_CurveClipboard, SIGNAL(clicked()), this, SLOT(ClipboardCurveButtonClicked()) ); + connect( m_Controls->m_ValuesClipboard, SIGNAL(clicked()), this, SLOT(ClipboardStatisticsButtonClicked()) ); + + } + + QString dstar = QString::number(m_Controls->m_DStarSlider->value()/1000.0); + m_Controls->m_DStarLabel->setText(dstar); + + QString bthresh = QString::number(m_Controls->m_BThreshSlider->value()*5.0); + m_Controls->m_BThreshLabel->setText(bthresh); + + QString s0thresh = QString::number(m_Controls->m_S0ThreshSlider->value()*0.5); + m_Controls->m_S0ThreshLabel->setText(s0thresh); + + QString numits = QString::number(m_Controls->m_NumItSlider->value()); + m_Controls->m_NumItsLabel->setText(numits); + + QString lambda = QString::number(m_Controls->m_LambdaSlider->value()*.00001); + m_Controls->m_LambdaLabel->setText(lambda); + + m_Controls->m_VisualizeResultsWidget->setVisible(m_Controls->m_DisplayResultsCheckbox->isChecked()); + m_Controls->m_MethodCombo->setVisible(m_Controls->m_ChooseMethod->isChecked()); + +// m_Controls->m_ADCBValues->setVisible(m_Controls->m_CheckADC->isChecked()); + + MethodCombo(m_Controls->m_MethodCombo->currentIndex()); + +} + +void QmitkIVIMView::Checkbox() +{ + m_Controls->m_VisualizeResultsWidget->setVisible(m_Controls->m_DisplayResultsCheckbox->isChecked()); + +// m_Controls->m_ADCBValues->setVisible(m_Controls->m_CheckADC->isChecked()); + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::MethodCombo(int val) +{ + switch(val) + { + case 0: + m_Controls->dstar->setVisible(false); + m_Controls->thres->setVisible(false); + m_Controls->thres_2->setVisible(true); + m_Controls->m_RegFrame->setVisible(false); + break; + case 1: + m_Controls->dstar->setVisible(true); + m_Controls->thres->setVisible(false); + m_Controls->thres_2->setVisible(true); + m_Controls->m_RegFrame->setVisible(false); + break; + case 2: + m_Controls->dstar->setVisible(false); + m_Controls->thres->setVisible(true); + m_Controls->thres_2->setVisible(true); + m_Controls->m_RegFrame->setVisible(false); + break; + case 3: + m_Controls->dstar->setVisible(false); + m_Controls->thres->setVisible(true); + m_Controls->thres_2->setVisible(true); + m_Controls->m_RegFrame->setVisible(false); + break; + case 4: + m_Controls->dstar->setVisible(false); + m_Controls->thres->setVisible(true); + m_Controls->thres_2->setVisible(true); + m_Controls->m_RegFrame->setVisible(true); + break; + } + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::DStarSlider (int val) +{ + QString sval = QString::number(val/1000.0); + m_Controls->m_DStarLabel->setText(sval); + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::BThreshSlider (int val) +{ + QString sval = QString::number(val*5.0); + m_Controls->m_BThreshLabel->setText(sval); + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::S0ThreshSlider (int val) +{ + QString sval = QString::number(val*0.5); + m_Controls->m_S0ThreshLabel->setText(sval); + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::NumItsSlider (int val) +{ + QString sval = QString::number(val); + m_Controls->m_NumItsLabel->setText(sval); + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::LambdaSlider (int val) +{ + QString sval = QString::number(val*.00001); + m_Controls->m_LambdaLabel->setText(sval); + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) +{ + m_MultiWidget = &stdMultiWidget; + + { + mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget1->GetSliceNavigationController(); + itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); + command->SetCallbackFunction( this, &QmitkIVIMView::OnSliceChanged ); + m_SliceObserverTag1 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(NULL, 0), command ); + } + + { + mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget2->GetSliceNavigationController(); + itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); + command->SetCallbackFunction( this, &QmitkIVIMView::OnSliceChanged ); + m_SliceObserverTag2 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(NULL, 0), command ); + } + + { + mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget3->GetSliceNavigationController(); + itk::ReceptorMemberCommand::Pointer command = itk::ReceptorMemberCommand::New(); + command->SetCallbackFunction( this, &QmitkIVIMView::OnSliceChanged ); + m_SliceObserverTag3 = slicer->AddObserver( mitk::SliceNavigationController::GeometrySliceEvent(NULL, 0), command ); + } + +} + +void QmitkIVIMView::StdMultiWidgetNotAvailable() +{ + + { + mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget1->GetSliceNavigationController(); + slicer->RemoveObserver( m_SliceObserverTag1 ); + } + + { + mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget2->GetSliceNavigationController(); + slicer->RemoveObserver( m_SliceObserverTag2 ); + } + + { + mitk::SliceNavigationController* slicer = m_MultiWidget->mitkWidget3->GetSliceNavigationController(); + slicer->RemoveObserver( m_SliceObserverTag3 ); + } + + m_MultiWidget = NULL; +} + +void QmitkIVIMView::OnSelectionChanged( std::vector nodes ) +{ + bool foundOneDiffusionImage = false; + + // iterate all selected objects, adjust warning visibility + for( std::vector::iterator it = nodes.begin(); + it != nodes.end(); + ++it ) + { + mitk::DataNode::Pointer node = *it; + + if( node.IsNotNull() ) + { + mitk::DiffusionImage* img = dynamic_cast*>(node->GetData()); + if( img ) + { + if(!foundOneDiffusionImage ) + { + foundOneDiffusionImage = true; + } + else + { + foundOneDiffusionImage = false; + } + } + } + } + +// m_Controls->m_ADCBValues->setVisible( foundOneDiffusionImage && m_Controls->m_CheckADC->isChecked() ); + + m_Controls->m_ButtonStart->setEnabled( foundOneDiffusionImage ); + m_Controls->m_ButtonAutoThres->setEnabled( foundOneDiffusionImage ); + + m_Controls->m_ControlsFrame->setEnabled( foundOneDiffusionImage ); + m_Controls->m_BottomControlsFrame->setEnabled( foundOneDiffusionImage ); + + itk::StartEvent dummy; + OnSliceChanged(dummy); +} + +void QmitkIVIMView::AutoThreshold() +{ + std::vector nodes = this->GetDataManagerSelection(); + if (nodes.empty()) return; + + if (!nodes.front()) + { + // Nothing selected. Inform the user and return + QMessageBox::information( NULL, "Template", "Please load and select a diffusion image before starting image processing."); + return; + } + + typedef mitk::DiffusionImage DiffImgType; + DiffImgType* dimg = dynamic_cast(nodes.front()->GetData()); + + if (!dimg) + { + // Nothing selected. Inform the user and return + QMessageBox::information( NULL, "Template", "No valid diffusion image was found."); + return; + } + + // find bzero index + int index = -1; + DiffImgType::GradientDirectionContainerType::Pointer directions = dimg->GetDirections(); + for(DiffImgType::GradientDirectionContainerType::ConstIterator it = directions->Begin(); + it != directions->End(); ++it) + { + index++; + DiffImgType::GradientDirectionType g = it.Value(); + if(g[0] == 0 && g[1] == 0 && g[2] == 0 ) + break; + } + + typedef itk::VectorImage VecImgType; + VecImgType::Pointer vecimg = dimg->GetVectorImage(); + + int vecLength = vecimg->GetVectorLength(); + index = index > vecLength-1 ? vecLength-1 : index; + + MITK_INFO << "Performing Histogram Analysis on Channel" << index; + + typedef itk::Image ImgType; + ImgType::Pointer img = ImgType::New(); + mitk::CastToItkImage(dimg, img); + + itk::ImageRegionIterator itw (img, img->GetLargestPossibleRegion() ); + itw = itw.Begin(); + + itk::ImageRegionConstIterator itr (vecimg, vecimg->GetLargestPossibleRegion() ); + itr = itr.Begin(); + + while(!itr.IsAtEnd()) + { + itw.Set(itr.Get().GetElement(index)); + ++itr; + ++itw; + } + + typedef itk::Statistics::ScalarImageToHistogramGenerator< ImgType > + HistogramGeneratorType; + typedef HistogramGeneratorType::HistogramType HistogramType; + + HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); + histogramGenerator->SetInput( img ); + histogramGenerator->SetMarginalScale( 10 ); // Defines y-margin width of histogram + histogramGenerator->SetNumberOfBins( 100 ); // CT range [-1024, +2048] --> bin size 4 values + histogramGenerator->SetHistogramMin( dimg->GetScalarValueMin() ); + histogramGenerator->SetHistogramMax( dimg->GetScalarValueMax() * .5 ); + histogramGenerator->Compute(); + + HistogramType::ConstIterator iter = histogramGenerator->GetOutput()->Begin(); + float maxFreq = 0; + float maxValue = 0; + while ( iter != histogramGenerator->GetOutput()->End() ) + { + if(iter.GetFrequency() > maxFreq) + { + maxFreq = iter.GetFrequency(); + maxValue = iter.GetMeasurementVector()[0]; + } + ++iter; + } + + maxValue *= 2; + + int sliderPos = maxValue * 2; + m_Controls->m_S0ThreshSlider->setValue(sliderPos); + S0ThreshSlider(sliderPos); +} + +void QmitkIVIMView::FittIVIMStart() +{ + + std::vector nodes = this->GetDataManagerSelection(); + if (nodes.empty()) return; + + if (!nodes.front()) + { + // Nothing selected. Inform the user and return + QMessageBox::information( NULL, "Template", "Please load and select a diffusion image before starting image processing."); + return; + } + + mitk::DiffusionImage* img = + dynamic_cast*>( + nodes.front()->GetData()); + + if (!img) + { + // Nothing selected. Inform the user and return + QMessageBox::information( NULL, "Template", "No valid diffusion image was found."); + return; + } + + typedef itk::VectorImage VecImgType; + VecImgType::Pointer vecimg = img->GetVectorImage(); + + OutImgType::IndexType dummy; + FittIVIM(vecimg, img->GetDirections(), img->GetB_Value(), true, dummy); + OutputToDatastorage(nodes); +} + +void QmitkIVIMView::OnSliceChanged(const itk::EventObject& /*e*/) +{ + m_Controls->m_VisualizeResultsWidget->setVisible(false); + m_Controls->m_Warning->setVisible(false); + + if(!m_Controls->m_DisplayResultsCheckbox->isChecked()) return; + + std::vector nodes = this->GetDataManagerSelection(); + if (nodes.empty()) return; + if (!nodes.front()) return; + if (nodes.size()>2) return; + + mitk::DiffusionImage* diffusionImg = 0; + mitk::DiffusionImage* img1 = + dynamic_cast*>( + nodes.front()->GetData()); + + mitk::DiffusionImage* img2 = 0; + mitk::Image* maskImg = 0; + if(nodes.size()>1) + { + if(img1) + { + if(strcmp(nodes.at(1)->GetData()->GetNameOfClass(), "Image") != 0 ) + return; + + maskImg = dynamic_cast( + nodes.at(1)->GetData()); + + diffusionImg = img1; + } + else + { + if(strcmp(nodes.front()->GetData()->GetNameOfClass(), "Image") != 0 ) + return; + + maskImg = dynamic_cast( + nodes.front()->GetData()); + + diffusionImg = dynamic_cast*>( + nodes.at(1)->GetData()); + } + } + else + { + diffusionImg = img1; + } + + if (nodes.size()==2 && (!diffusionImg || !maskImg || m_Controls->m_MethodCombo->currentIndex() == 4 )) return; + if (nodes.size()==1 && !diffusionImg) return; + + if (!m_MultiWidget) return; + + m_Controls->m_VisualizeResultsWidget->setVisible(true); + + typedef itk::VectorImage VecImgType; + VecImgType::Pointer vecimg = (VecImgType*)diffusionImg->GetVectorImage().GetPointer(); + + VecImgType::Pointer roiImage = VecImgType::New(); + if(maskImg == 0) + { + int roisize = 0; + if(m_Controls->m_MethodCombo->currentIndex() == 4) + roisize = 5; + + mitk::Point3D pos = m_MultiWidget->GetCrossPosition(); + VecImgType::IndexType crosspos; + diffusionImg->GetTimeSlicedGeometry()->WorldToIndex(pos, crosspos); + + VecImgType::IndexType index; + index[0] = crosspos[0] - roisize; index[0] = index[0] < 0 ? 0 : index[0]; + index[1] = crosspos[1] - roisize; index[1] = index[1] < 0 ? 0 : index[1]; + index[2] = crosspos[2] - roisize; index[2] = index[2] < 0 ? 0 : index[2]; + + VecImgType::SizeType size; + size[0] = roisize*2+1; + size[1] = roisize*2+1; + size[2] = roisize*2+1; + + VecImgType::SizeType maxSize = vecimg->GetLargestPossibleRegion().GetSize(); + size[0] = index[0]+size[0] > maxSize[0] ? maxSize[0]-index[0] : size[0]; + size[1] = index[1]+size[1] > maxSize[1] ? maxSize[1]-index[1] : size[1]; + size[2] = index[2]+size[2] > maxSize[2] ? maxSize[2]-index[2] : size[2]; + + VecImgType::RegionType region; + region.SetSize( size ); + region.SetIndex( index ); + vecimg->SetRequestedRegion( region ); + + VecImgType::IndexType newstart; + newstart.Fill(0); + + VecImgType::RegionType newregion; + newregion.SetSize( size ); + newregion.SetIndex( newstart ); + + roiImage->CopyInformation( vecimg ); + roiImage->SetRegions( newregion ); + roiImage->SetOrigin( pos ); + roiImage->Allocate(); + roiImage->SetPixel(newstart, vecimg->GetPixel(index)); + + FittIVIM(roiImage, diffusionImg->GetDirections(), diffusionImg->GetB_Value(), false, crosspos); + } + else + { + typedef itk::Image MaskImgType; + + MaskImgType::Pointer maskItk; + CastToItkImage( maskImg, maskItk ); + + mitk::Point3D pos; + pos[0] = 0; + pos[1] = 0; + pos[2] = 0; + + VecImgType::IndexType index; + index[0] = 0; + index[1] = 0; + index[2] = 0; + + VecImgType::SizeType size; + size[0] = 1; + size[1] = 1; + size[2] = 1; + + VecImgType::RegionType region; + region.SetSize( size ); + region.SetIndex( index ); + vecimg->SetRequestedRegion( region ); + + // iterators over output and input + itk::ImageRegionConstIteratorWithIndex + vecit(vecimg, vecimg->GetLargestPossibleRegion()); + + itk::VariableLengthVector avg(vecimg->GetVectorLength()); + avg.Fill(0); + + float numPixels = 0; + while ( ! vecit.IsAtEnd() ) + { + VecImgType::PointType point; + vecimg->TransformIndexToPhysicalPoint(vecit.GetIndex(), point); + + MaskImgType::IndexType index; + maskItk->TransformPhysicalPointToIndex(point, index); + + if(maskItk->GetPixel(index) != 0) + { + avg += vecit.Get(); + numPixels += 1.0; + } + + // update iterators + ++vecit; + + } + + avg /= numPixels; + + m_Controls->m_Warning->setText(QString("Averaging ")+QString::number((int)numPixels)+QString(" voxels!")); + m_Controls->m_Warning->setVisible(true); + + roiImage->CopyInformation( vecimg ); + roiImage->SetRegions( region ); + roiImage->SetOrigin( pos ); + roiImage->Allocate(); + roiImage->SetPixel(index, avg); + + FittIVIM(roiImage, diffusionImg->GetDirections(), diffusionImg->GetB_Value(), false, index); + } + + vecimg->SetRegions( vecimg->GetLargestPossibleRegion() ); + + m_Controls->m_VisualizeResultsWidget->SetParameters(m_Snap); + +} + +void QmitkIVIMView::FittIVIM(itk::VectorImage* vecimg, DirContainerType* dirs, float bval, bool multivoxel, OutImgType::IndexType &crosspos) +{ + + IVIMFilterType::Pointer filter = IVIMFilterType::New(); + filter->SetInput(vecimg); + filter->SetGradientDirections(dirs); + filter->SetBValue(bval); + + switch(m_Controls->m_MethodCombo->currentIndex()) + { + + case 0: + filter->SetMethod(IVIMFilterType::IVIM_FIT_ALL); + filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); + break; + + case 1: + filter->SetMethod(IVIMFilterType::IVIM_DSTAR_FIX); + filter->SetDStar(m_Controls->m_DStarLabel->text().toDouble()); + filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); + break; + + case 2: + filter->SetMethod(IVIMFilterType::IVIM_D_THEN_DSTAR); + filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); + filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); + filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); + break; + + case 3: + filter->SetMethod(IVIMFilterType::IVIM_LINEAR_D_THEN_F); + filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); + filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); + filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); + break; + + case 4: + filter->SetMethod(IVIMFilterType::IVIM_REGULARIZED); + filter->SetBThres(m_Controls->m_BThreshLabel->text().toDouble()); + filter->SetS0Thres(m_Controls->m_S0ThreshLabel->text().toDouble()); + filter->SetNumberIterations(m_Controls->m_NumItsLabel->text().toInt()); + filter->SetLambda(m_Controls->m_LambdaLabel->text().toDouble()); + filter->SetFitDStar(m_Controls->m_CheckDStar->isChecked()); + break; + } + + if(!multivoxel) + { + filter->SetFitDStar(true); + } + + filter->SetNumberOfThreads(1); + filter->SetVerbose(multivoxel); + filter->SetCrossPosition(crosspos); + filter->Update(); + + m_Snap = filter->GetSnapshot(); + m_DStarMap = filter->GetOutput(2); + m_DMap = filter->GetOutput(1); + m_fMap = filter->GetOutput(0); + +} + +void QmitkIVIMView::OutputToDatastorage(std::vector nodes) +{ + // Outputs to Datastorage + QString basename(nodes.front()->GetName().c_str()); + + if(m_Controls->m_CheckDStar->isChecked()) + { + mitk::Image::Pointer dstarimage = mitk::Image::New(); + dstarimage->InitializeByItk(m_DStarMap.GetPointer()); + dstarimage->SetVolume(m_DStarMap->GetBufferPointer()); + QString newname2 = basename; newname2 = newname2.append("_DStarMap%1").arg(m_Controls->m_MethodCombo->currentIndex()); + mitk::DataNode::Pointer node2=mitk::DataNode::New(); + node2->SetData( dstarimage ); + node2->SetName(newname2.toAscii()); + GetDefaultDataStorage()->Add(node2); + } + + if(m_Controls->m_CheckD->isChecked()) + { + mitk::Image::Pointer dimage = mitk::Image::New(); + dimage->InitializeByItk(m_DMap.GetPointer()); + dimage->SetVolume(m_DMap->GetBufferPointer()); + QString newname1 = basename; newname1 = newname1.append("_DMap%1").arg(m_Controls->m_MethodCombo->currentIndex()); + mitk::DataNode::Pointer node1=mitk::DataNode::New(); + node1->SetData( dimage ); + node1->SetName(newname1.toAscii()); + GetDefaultDataStorage()->Add(node1); + } + + if(m_Controls->m_Checkf->isChecked()) + { + mitk::Image::Pointer image = mitk::Image::New(); + image->InitializeByItk(m_fMap.GetPointer()); + image->SetVolume(m_fMap->GetBufferPointer()); + QString newname0 = basename; newname0 = newname0.append("_fMap%1").arg(m_Controls->m_MethodCombo->currentIndex()); + mitk::DataNode::Pointer node=mitk::DataNode::New(); + node->SetData( image ); + node->SetName(newname0.toAscii()); + GetDefaultDataStorage()->Add(node); + } + + m_MultiWidget->RequestUpdate(); + +} + +void QmitkIVIMView::ChooseMethod() +{ + m_Controls->m_MethodCombo->setVisible(m_Controls->m_ChooseMethod->isChecked()); +} + +void QmitkIVIMView::ClipboardCurveButtonClicked() +{ + if(true) + { + + QString clipboard("Measurement Points\n"); + for ( int i=0; isetText( + clipboard, QClipboard::Clipboard ); + } + else + { + QApplication::clipboard()->clear(); + } +} + + +void QmitkIVIMView::ClipboardStatisticsButtonClicked() +{ + if ( true ) + { + + QString clipboard( "f \t D \t D* \n" ); + clipboard = clipboard.append( "%L1 \t %L2 \t %L3" ) + .arg( m_Snap.currentF, 0, 'f', 10 ) + .arg( m_Snap.currentD, 0, 'f', 10 ) + .arg( m_Snap.currentDStar, 0, 'f', 10 ) ; + + QApplication::clipboard()->setText( + clipboard, QClipboard::Clipboard ); + } + else + { + QApplication::clipboard()->clear(); + } +} diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.h new file mode 100644 index 0000000000..8ec0ab6db1 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMView.h @@ -0,0 +1,108 @@ +/*========================================================================= +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 _QMITKIVIMVIEW_H_INCLUDED +#define _QMITKIVIMVIEW_H_INCLUDED + +#include + +#include + +#include "ui_QmitkIVIMViewControls.h" + +#include "itkVectorImage.h" +#include "itkImage.h" + +#include "mitkDiffusionImage.h" +#include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" + +/*! + \brief QmitkIVIMView + + \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 +*/ + +class QmitkIVIMView : 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: + + static const std::string VIEW_ID; + + QmitkIVIMView(); + virtual ~QmitkIVIMView(); + + typedef mitk::DiffusionImage::GradientDirectionContainerType DirContainerType; + typedef itk::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter IVIMFilterType; + typedef itk::Image OutImgType; + + virtual void CreateQtPartControl(QWidget *parent); + + virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); + virtual void StdMultiWidgetNotAvailable(); + + void OnSliceChanged(const itk::EventObject& e); + void OutputToDatastorage(std::vector nodes); + void FittIVIM(itk::VectorImage* vecimg, DirContainerType* dirs, float bval, bool multivoxel, OutImgType::IndexType &crosspos); + +protected slots: + + /// \brief Called when the user clicks the GUI button + void FittIVIMStart(); + void AutoThreshold(); + + void MethodCombo(int val); + void Checkbox(); + + void DStarSlider(int val); + void BThreshSlider(int val); + void S0ThreshSlider(int val); + void NumItsSlider(int val); + void LambdaSlider(int val); + void ChooseMethod(); + void ClipboardStatisticsButtonClicked(); + void ClipboardCurveButtonClicked(); + +protected: + + /// \brief called by QmitkFunctionality when DataManager's selection has changed + virtual void OnSelectionChanged( std::vector nodes ); + + Ui::QmitkIVIMViewControls* m_Controls; + + QmitkStdMultiWidget* m_MultiWidget; + int m_SliceObserverTag1; + int m_SliceObserverTag2; + int m_SliceObserverTag3; + + OutImgType::Pointer m_DStarMap; + OutImgType::Pointer m_DMap; + OutImgType::Pointer m_fMap; + + IVIMFilterType::IVIMSnapshot m_Snap; + +}; + + + +#endif // _QMITKIVIMVIEW_H_INCLUDED + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMViewControls.ui new file mode 100644 index 0000000000..baa5bbde7c --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkIVIMViewControls.ui @@ -0,0 +1,613 @@ + + + QmitkIVIMViewControls + + + + 0 + 0 + 307 + 790 + + + + + 0 + 0 + + + + QmitkTemplate + + + + 0 + + + + + Intra Voxel Incoherent Motion Estimation + + + + 0 + + + 9 + + + 0 + + + 0 + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 80 + 16777215 + + + + D* + + + + + + + 100 + + + 60 + + + Qt::Horizontal + + + + + + + + 51 + 16777215 + + + + 200 + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 80 + 16777215 + + + + neglect b< + + + + + + + 250 + + + 34 + + + Qt::Horizontal + + + + + + + + 51 + 16777215 + + + + 46.5 + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 80 + 16777215 + + + + neglect Si< + + + + + + + 100 + + + 0 + + + Qt::Horizontal + + + + + + + + 30 + 16777215 + + + + TextLabel + + + + + + + + 15 + 16777215 + + + + * + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 80 + 16777215 + + + + #iterations + + + + + + + 100 + + + 10 + + + Qt::Horizontal + + + + + + + + 30 + 16777215 + + + + TextLabel + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 80 + 16777215 + + + + lambda + + + + + + + 1000 + + + 10 + + + Qt::Horizontal + + + + + + + + 30 + 16777215 + + + + TextLabel + + + + + + + + 15 + 16777215 + + + + * + + + + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 80 + 0 + + + + Output Images + + + + + + + f + + + true + + + + + + + D + + + false + + + + + + + D* + + + false + + + + + + + + + + Generate Output Images + + + + + + + color: rgb(255, 0, 0); +font: 75 14pt "Ubuntu"; + + + Bla bla bla + + + Qt::RichText + + + + + + + display voxel-wise results + + + true + + + + + + + + + + true + + + + 0 + 0 + + + + + 0 + 400 + + + + + 0 + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + Curve to Clipboard + + + + + + + Values to Clipboard + + + + + + + + + + Choose Method + + + + + + + 2 + + + + 3 Param. Fit + + + + + Fit D & f with fixed D* value + + + + + Fit D & f (high b), then fit D* + + + + + Linearly fit D & f (high b), then fit D* + + + + + Regularized + + + + + + + + + + + + + + Qt::Vertical + + + QSizePolicy::Expanding + + + + 20 + 220 + + + + + + + + + + QmitkIVIMWidget + QWidget +
QmitkIVIMWidget.h
+ 1 + + + + + diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp new file mode 100644 index 0000000000..6cc92bd4fa --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp @@ -0,0 +1,768 @@ +/*======================================================================= +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 +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp +#define __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp + +#include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h" +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionIterator.h" + +#include "vnl/vnl_matrix.h" +#include "vnl/algo/vnl_symmetric_eigensystem.h" + +#include "itkRegularizedIVIMReconstructionFilter.h" + +#include + +#define IVIM_FOO -100000 + +namespace itk { + + + template< class TIn, class TOut> + DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() : + m_GradientDirectionContainer(NULL), + m_Method(IVIM_DSTAR_FIX), + m_FitDStar(true), + m_Verbose(true) + { + this->SetNumberOfRequiredInputs( 1 ); + + this->SetNumberOfRequiredOutputs( 3 ); + typename OutputImageType::Pointer outputPtr1 = OutputImageType::New(); + this->SetNthOutput(0, outputPtr1.GetPointer()); + typename OutputImageType::Pointer outputPtr2 = OutputImageType::New(); + this->SetNthOutput(1, outputPtr2.GetPointer()); + typename OutputImageType::Pointer outputPtr3 = OutputImageType::New(); + this->SetNthOutput(2, outputPtr3.GetPointer()); + } + + + + template< class TIn, class TOut> + void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::BeforeThreadedGenerateData() + { + + // Input must be an itk::VectorImage. + std::string gradientImageClassName( + this->ProcessObject::GetInput(0)->GetNameOfClass()); + if ( strcmp(gradientImageClassName.c_str(),"VectorImage") != 0 ) + { + itkExceptionMacro( << + "There is only one Gradient image. I expect that to be a VectorImage. " + << "But its of type: " << gradientImageClassName ); + } + + typename InputImageType::Pointer img = static_cast< InputImageType * >( + this->ProcessObject::GetInput(0) ); + + // Compute the indicies of the baseline images and gradient images + GradientDirectionContainerType::ConstIterator gdcit = + this->m_GradientDirectionContainer->Begin(); + while( gdcit != this->m_GradientDirectionContainer->End() ) + { + if(gdcit.Value().one_norm() <= 0.0) + { + m_Snap.baselineind.push_back(gdcit.Index()); + } + else + { + m_Snap.gradientind.push_back(gdcit.Index()); + double twonorm = gdcit.Value().two_norm(); + m_Snap.bvals.push_back( m_BValue*twonorm*twonorm ); + } + + ++gdcit; + } + + + // check sind die grad und base gleichlang? alle grad gerade und base ungerade? dann iterierende aufnahme!! + m_Snap.iterated_sequence = false; + if(m_Snap.baselineind.size() == m_Snap.gradientind.size()) + { + int size = m_Snap.baselineind.size(); + int sum_b = 0, sum_g = 0; + for(int i=0; im_BThres) + { + m_Snap.high_indices.push_back(i); + } + } + } + m_Snap.Nhigh = m_Snap.high_indices.size(); + m_Snap.high_bvalues.set_size(m_Snap.Nhigh); + m_Snap.high_meas.set_size(m_Snap.Nhigh); + for(int i=0; i + MeasAndBvals DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::ApplyS0Threshold(vnl_vector &meas, vnl_vector &bvals) + { + std::vector newmeas; + std::vector newbvals; + + int N = meas.size(); + for(int i=0; i + void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int ) + { + + typename OutputImageType::Pointer outputImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); + oit.GoToBegin(); + + typename OutputImageType::Pointer dImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); + ImageRegionIterator< OutputImageType > oit1(dImage, outputRegionForThread); + oit1.GoToBegin(); + + typename OutputImageType::Pointer dstarImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2)); + ImageRegionIterator< OutputImageType > oit2(dstarImage, outputRegionForThread); + oit2.GoToBegin(); + + typedef ImageRegionConstIterator< InputImageType > InputIteratorType; + typedef typename InputImageType::PixelType InputVectorType; + typename InputImageType::Pointer inputImagePointer = NULL; + + // Would have liked a dynamic_cast here, but seems SGI doesn't like it + // The enum will DiffusionIntravoxelIncoherentMotionReconstructionImageFilterensure that an inappropriate cast is not done + inputImagePointer = static_cast< InputImageType * >( + this->ProcessObject::GetInput(0) ); + + InputIteratorType iit(inputImagePointer, outputRegionForThread ); + iit.GoToBegin(); + + // init internal vector image for regularized fit + m_InternalVectorImage = VectorImageType::New(); + m_InternalVectorImage->SetSpacing( inputImagePointer->GetSpacing() ); // Set the image spacing + m_InternalVectorImage->SetOrigin( inputImagePointer->GetOrigin() ); // Set the image origin + m_InternalVectorImage->SetDirection( inputImagePointer->GetDirection() ); // Set the image direction + m_InternalVectorImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() ); + + m_InitialFitImage = InitialFitImageType::New(); + m_InitialFitImage->SetSpacing( inputImagePointer->GetSpacing() ); // Set the image spacing + m_InitialFitImage->SetOrigin( inputImagePointer->GetOrigin() ); // Set the image origin + m_InitialFitImage->SetDirection( inputImagePointer->GetDirection() ); // Set the image direction + m_InitialFitImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() ); + + if(m_Method == IVIM_REGULARIZED) + { + m_InternalVectorImage->SetVectorLength(m_Snap.Nhigh); + m_InternalVectorImage->Allocate(); + VectorImageType::PixelType varvec(m_Snap.Nhigh); + for(int i=0; iFillBuffer(varvec); + + m_InitialFitImage->Allocate(); + InitialFitImageType::PixelType vec; + vec[0] = 0.5; vec[1] = 0.01; vec[2]=0.001; + m_InitialFitImage->FillBuffer(vec); + } + + typedef itk::ImageRegionIterator VectorIteratorType; + VectorIteratorType vecit(m_InternalVectorImage, outputRegionForThread ); + vecit.GoToBegin(); + + typedef itk::ImageRegionIterator InitIteratorType; + InitIteratorType initit(m_InitialFitImage, outputRegionForThread ); + initit.GoToBegin(); + + while( !iit.IsAtEnd() ) + { + InputVectorType measvec = iit.Get(); + + typename NumericTraits::AccumulateType b0 = NumericTraits::Zero; + + m_Snap.meas.set_size(m_Snap.N); + m_Snap.allmeas.set_size(m_Snap.N); + if(!m_Snap.iterated_sequence) + { + // Average the baseline image pixels + for(unsigned int i = 0; i < m_Snap.baselineind.size(); ++i) + { + b0 += measvec[m_Snap.baselineind[i]]; + } + b0 /= m_Snap.baselineind.size(); + + // measurement vector + for(int i = 0; i < m_Snap.N; ++i) + { + m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); + + if(measvec[m_Snap.gradientind[i]] > m_S0Thres) + { + m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); + } + else + { + m_Snap.meas[i] = IVIM_FOO; + } + } + } + else + { + // measurement vector + for(int i = 0; i < m_Snap.N; ++i) + { + b0 = measvec[m_Snap.baselineind[i]]; + + m_Snap.allmeas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); + + if(measvec[m_Snap.gradientind[i]] > m_S0Thres) + { + m_Snap.meas[i] = measvec[m_Snap.gradientind[i]] / (b0+.0001); + } + else + { + m_Snap.meas[i] = IVIM_FOO; + } + } + } + + m_Snap.currentF = 0; + m_Snap.currentD = 0; + m_Snap.currentDStar = 0; + + switch(m_Method) + { + + case IVIM_D_THEN_DSTAR: + { + + for(int i=0; i x_donly(2); + x_donly[0] = 0.001; + x_donly[1] = 0.1; + // f 0.1 Dstar 0.01 D 0.001 + + vnl_levenberg_marquardt lm_donly(f_donly); + lm_donly.set_f_tolerance(0.0001); + lm_donly.minimize(x_donly); + m_Snap.currentD = x_donly[0]; + m_Snap.currentF = x_donly[1]; + + + if(m_FitDStar) + { + MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + m_Snap.bvals2 = input2.bvals; + m_Snap.meas2 = input2.meas; + if (input2.N < 2) break; + + IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF); + f_dstar_only.set_bvalues(input2.bvals); + f_dstar_only.set_measurements(input2.meas); + + vnl_vector< double > x_dstar_only(1); + vnl_vector< double > fx_dstar_only(input2.N); + + double opt = 1111111111111111.0; + int opt_idx = -1; + int num_its = 100; + double min_val = .001; + double max_val = .15; + for(int i=0; i x_fixd(2); + // x_fixd[0] = 0.1; + // x_fixd[1] = 0.01; + // // f 0.1 Dstar 0.01 D 0.001 + + // vnl_levenberg_marquardt lm_fixd(f_fixd); + // lm_fixd.set_f_tolerance(0.0001); + // lm_fixd.minimize(x_fixd); + + // m_Snap.currentF = x_fixd[0]; + // m_Snap.currentDStar = x_fixd[1]; + } + + break; + } + + case IVIM_DSTAR_FIX: + { + MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + m_Snap.bvals1 = input.bvals; + m_Snap.meas1 = input.meas; + if (input.N < 2) break; + + IVIM_fixdstar f_fixdstar(input.N,m_DStar); + f_fixdstar.set_bvalues(input.bvals); + f_fixdstar.set_measurements(input.meas); + + vnl_vector< double > x(2); + x[0] = 0.1; + x[1] = 0.001; + // f 0.1 Dstar 0.01 D 0.001 + + vnl_levenberg_marquardt lm(f_fixdstar); + lm.set_f_tolerance(0.0001); + lm.minimize(x); + + m_Snap.currentF = x[0]; + m_Snap.currentD = x[1]; + m_Snap.currentDStar = m_DStar; + + break; + } + + case IVIM_FIT_ALL: + { + + MeasAndBvals input = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + m_Snap.bvals1 = input.bvals; + m_Snap.meas1 = input.meas; + if (input.N < 3) break; + + IVIM_3param f_3param(input.N); + f_3param.set_bvalues(input.bvals); + f_3param.set_measurements(input.meas); + + vnl_vector< double > x(3); + x[0] = 0.1; + x[1] = 0.001; + x[2] = 0.01; + // f 0.1 Dstar 0.01 D 0.001 + + vnl_levenberg_marquardt lm(f_3param); + lm.set_f_tolerance(0.0001); + lm.minimize(x); + + m_Snap.currentF = x[0]; + m_Snap.currentD = x[1]; + m_Snap.currentDStar = x[2]; + + break; + } + + case IVIM_LINEAR_D_THEN_F: + { + + // // neglect zero-measurements + // bool zero = false; + // for(int i=0; i X(input.N,2); + for(int i=0; i XX = X.transpose() * X; + vnl_symmetric_eigensystem eigs(XX); + + vnl_vector eig; + if(eigs.get_eigenvalue(0) > eigs.get_eigenvalue(1)) + eig = eigs.get_eigenvector(0); + else + eig = eigs.get_eigenvector(1); + + m_Snap.currentF = 1 - exp( meas_m - bval_m*(eig(1)/eig(0)) ); + m_Snap.currentD = -eig(1)/eig(0); + + if(m_FitDStar) + { + MeasAndBvals input2 = ApplyS0Threshold(m_Snap.meas, m_Snap.bvalues); + m_Snap.bvals2 = input2.bvals; + m_Snap.meas2 = input2.meas; + if (input2.N < 2) break; + + IVIM_dstar_only f_dstar_only(input2.N,m_Snap.currentD,m_Snap.currentF); + f_dstar_only.set_bvalues(input2.bvals); + f_dstar_only.set_measurements(input2.meas); + + vnl_vector< double > x_dstar_only(1); + vnl_vector< double > fx_dstar_only(input2.N); + + double opt = 1111111111111111.0; + int opt_idx = -1; + int num_its = 100; + double min_val = .001; + double max_val = .15; + for(int i=0; i " << DStar; + // x_dstar_only[0] = 0.01; + // // f 0.1 Dstar 0.01 D 0.001 + + // vnl_levenberg_marquardt lm_dstar_only(f_dstar_only); + // lm_dstar_only.set_f_tolerance(0.0001); + // lm_dstar_only.minimize(x_dstar_only); + + // DStar = x_dstar_only[0]; + + break; + } + + case IVIM_REGULARIZED: + { + //m_Snap.high_meas, m_Snap.high_bvalues; + for(int i=0; i x_donly(2); + x_donly[0] = 0.001; + x_donly[1] = 0.1; + + if(input.N >= 2) + { + IVIM_d_and_f f_donly(input.N); + f_donly.set_bvalues(input.bvals); + f_donly.set_measurements(input.meas); + //MITK_INFO << "initial fit N=" << input.N << ", min-b = " << input.bvals[0] << ", max-b = " << input.bvals[input.N-1]; + vnl_levenberg_marquardt lm_donly(f_donly); + lm_donly.set_f_tolerance(0.0001); + lm_donly.minimize(x_donly); + } + + typename InitialFitImageType::PixelType initvec; + initvec[0] = x_donly[1]; + initvec[1] = x_donly[0]; + initit.Set(initvec); + //MITK_INFO << "Init vox " << initit.GetIndex() << " with " << initvec[0] << "; " << initvec[1]; + ++initit; + + int N = m_Snap.high_meas.size(); + typename VectorImageType::PixelType vec(N); + for(int i=0; i + void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::AfterThreadedGenerateData() + { + if(m_Method == IVIM_REGULARIZED) + { + typename OutputImageType::Pointer outputImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + ImageRegionIterator< OutputImageType > oit0(outputImage, outputImage->GetLargestPossibleRegion()); + oit0.GoToBegin(); + + typename OutputImageType::Pointer dImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(1)); + ImageRegionIterator< OutputImageType > oit1(dImage, dImage->GetLargestPossibleRegion()); + oit1.GoToBegin(); + + typename OutputImageType::Pointer dstarImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(2)); + ImageRegionIterator< OutputImageType > oit2(dstarImage, dstarImage->GetLargestPossibleRegion()); + oit2.GoToBegin(); + + typedef itk::RegularizedIVIMReconstructionFilter + RegFitType; + RegFitType::Pointer filter = RegFitType::New(); + filter->SetInput(m_InitialFitImage); + filter->SetReferenceImage(m_InternalVectorImage); + filter->SetBValues(m_Snap.high_bvalues); + filter->SetNumberIterations(m_NumberIterations); + filter->SetNumberOfThreads(1); + filter->SetLambda(m_Lambda); + filter->Update(); + typename RegFitType::OutputImageType::Pointer outimg = filter->GetOutput(); + + ImageRegionConstIterator< RegFitType::OutputImageType > iit(outimg, outimg->GetLargestPossibleRegion()); + iit.GoToBegin(); + + while( !iit.IsAtEnd() ) + { + double f = iit.Get()[0]; + IVIM_CEIL( f, 0.0, 1.0 ); + + oit0.Set( myround(f * 100.0) ); + oit1.Set( myround(iit.Get()[1] * 10000.0) ); + oit2.Set( myround(iit.Get()[2] * 1000.0) ); + + if(!m_Verbose) + { + // report the middle voxel + if( iit.GetIndex()[0] == m_CrossPosition[0] + && iit.GetIndex()[1] == m_CrossPosition[1] + && iit.GetIndex()[2] == m_CrossPosition[2] ) + { + m_Snap.currentF = f; + m_Snap.currentD = iit.Get()[1]; + m_Snap.currentDStar = iit.Get()[2]; + m_Snap.allmeas = m_tmp_allmeas; + MITK_INFO << "setting " << f << ";" << iit.Get()[1] << ";" << iit.Get()[2]; + } + } + + ++oit0; + ++oit1; + ++oit2; + ++iit; + } + } + } + + template< class TIn, class TOut> + double DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::myround(double number) + { + return number < 0.0 ? ceil(number - 0.5) : floor(number + 0.5); + } + + template< class TIn, class TOut> + void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::SetGradientDirections( GradientDirectionContainerType *gradientDirection ) + { + this->m_GradientDirectionContainer = gradientDirection; + this->m_NumberOfGradientDirections = gradientDirection->Size(); + } + + template< class TIn, class TOut> + void DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + ::PrintSelf(std::ostream& os, Indent indent) const + { + Superclass::PrintSelf(os,indent); + + if ( m_GradientDirectionContainer ) + { + os << indent << "GradientDirectionContainer: " + << m_GradientDirectionContainer << std::endl; + } + else + { + os << indent << + "GradientDirectionContainer: (Gradient directions not set)" << std::endl; + } + } + +} + +#endif // __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_cpp diff --git a/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h new file mode 100644 index 0000000000..02ee5e4bd3 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h @@ -0,0 +1,401 @@ +/*========================================================================= + + Program: Insight Segmentation & Registration Toolkit + Module: $RCSfile: itkDiffusionTensor3DReconstructionImageFilter.h,v $ + Language: C++ + Date: $Date: 2006-03-27 17:01:06 $ + Version: $Revision: 1.12 $ + + 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 + PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h +#define __itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h + +#include "itkImageToImageFilter.h" +//#include "vnl/vnl_matrix.h" +#include "vnl/vnl_vector_fixed.h" +#include "vnl/vnl_matrix.h" +#include "vnl/algo/vnl_svd.h" +#include "itkVectorContainer.h" +#include "itkVectorImage.h" +//#include "QuadProg.h" +#include "itkVectorImage.h" + +#include "vnl/vnl_least_squares_function.h" +#include "vnl/algo/vnl_levenberg_marquardt.h" +#include "vnl/vnl_math.h" + +#define IVIM_CEIL(val,u,o) (val) = \ + ( (val) < (u) ) ? ( (u) ) : ( ( (val)>(o) ) ? ( (o) ) : ( (val) ) ); + +namespace itk{ + + /////////////////////////////////////////////////////////////////////// + // baseclass for IVIM fitting algorithms + struct IVIM_base + { + + void set_measurements(const vnl_vector& x) + { + measurements.set_size(x.size()); + measurements.copy_in(x.data_block()); + } + + void set_bvalues(const vnl_vector& x) + { + bvalues.set_size(x.size()); + bvalues.copy_in(x.data_block()); + } + + vnl_vector measurements; + vnl_vector bvalues; + + int N; + + }; + + /////////////////////////////////////////////////////////////////////// + // fitt all three parameters + struct IVIM_3param : public IVIM_base, vnl_least_squares_function + { + + IVIM_3param(unsigned int number_of_measurements) : + vnl_least_squares_function(3, number_of_measurements, no_gradient) + { + N = get_number_of_residuals(); + } + + void f(const vnl_vector& x, vnl_vector& fx) { + + double ef = x[0]; + double D = x[1]; + double Dstar = x[2]; + + for(int s=0; s& x, vnl_vector& fx) { + + double ef = x[0]; + double D = x[1]; + + for(int s=0; s& x, vnl_vector& fx) { + + double D = x[0]; + double f = x[1]; + + for(int s=0; s& x, vnl_vector& fx) { + + double ef = x[0]; + double Dstar = x[1]; + + for(int s=0; s& x, vnl_vector& fx) { + + double Dstar = x[0]; + + for(int s=0; s meas; + vnl_vector bvals; + int N; + }; + + /** \class DiffusionIntravoxelIncoherentMotionReconstructionImageFilter + * + */ + + template< class TInputPixelType, + class TOutputPixelType> + class DiffusionIntravoxelIncoherentMotionReconstructionImageFilter : + public ImageToImageFilter< VectorImage< TInputPixelType, 3 >, + Image< TOutputPixelType, 3 > > + { + + public: + + struct IVIMSnapshot + { + double currentF; + double currentFunceiled; + double currentD; + double currentDStar; + + bool iterated_sequence; // wether each measurement has its own b0-acqu. + std::vector baselineind; // baseline image indicies + + int N; // total number of measurements + std::vector gradientind; // gradient image indicies + std::vector bvals; // b-values != 0 + vnl_vector bvalues; // copy of bvalues != 0 + vnl_vector meas; // all measurements, thresholded blanked out + vnl_vector allmeas; // all measurements + + int Nhigh; // number of used measurements + std::vector high_indices; // indices of used measurements + vnl_vector high_bvalues; // bvals of used measurements + vnl_vector high_meas; // used measurements + + vnl_vector meas1; + vnl_vector bvals1; + + vnl_vector meas2; + vnl_vector bvals2; + +// double currentADC; +// int Nadc; // number used measurements for ADC calculation +// std::vector adc_indices; // indices used for ADC calculation +// vnl_vector adc_bvalues; // bvals used for ADC calculation +// vnl_vector adc_meas; // for ADC calculation + + }; + + enum IVIM_Method + { + IVIM_FIT_ALL, + IVIM_DSTAR_FIX, + IVIM_D_THEN_DSTAR, + IVIM_LINEAR_D_THEN_F, + IVIM_REGULARIZED + }; + + typedef DiffusionIntravoxelIncoherentMotionReconstructionImageFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< VectorImage< TInputPixelType, 3>, + Image< TOutputPixelType,3 > > Superclass; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(DiffusionIntravoxelIncoherentMotionReconstructionImageFilter, + ImageToImageFilter); + + typedef TOutputPixelType OutputPixelType; + + typedef TInputPixelType InputPixelType; + + /** Reference image data, This image is aquired in the absence + * of a diffusion sensitizing field gradient */ + typedef typename Superclass::InputImageType InputImageType; + + typedef Image< OutputPixelType, 3 > OutputImageType; + + typedef typename Superclass::OutputImageRegionType + OutputImageRegionType; + + /** Holds each magnetic field gradient used to acquire one DWImage */ + typedef vnl_vector_fixed< double, 3 > GradientDirectionType; + + /** Container to hold gradient directions of the 'n' DW measurements */ + typedef VectorContainer< unsigned int, + GradientDirectionType > GradientDirectionContainerType; + + // vector image typedefs for regularized fit + typedef itk::VectorImage VectorImageType; + typedef itk::Image, 3> InitialFitImageType; + + /** set method to add gradient directions and its corresponding + * image. The image here is a VectorImage. The user is expected to pass the + * gradient directions in a container. The ith element of the container + * corresponds to the gradient direction of the ith component image the + * VectorImage. For the baseline image, a vector of all zeros + * should be set.*/ + void SetGradientDirections( GradientDirectionContainerType * ); + + void SetBValue(double bval){m_BValue = bval;} + void SetBThres(double bval){m_BThres = bval;} + void SetS0Thres(double val){m_S0Thres = val;} + void SetDStar(double dstar){m_DStar = dstar;} + void SetFitDStar(bool fit){m_FitDStar = fit;} + void SetVerbose(bool verbose){m_Verbose = verbose;} + void SetNumberIterations(int num){m_NumberIterations = num;} + void SetLambda(double lambda){m_Lambda = lambda;} + void SetCrossPosition(typename InputImageType::IndexType crosspos){this->m_CrossPosition = crosspos;} + void SetMethod(IVIM_Method method){m_Method = method;} +// void SetGradientIndicesForADCCalculation(std::vector inds){this->m_ADCInds = inds;} + + IVIMSnapshot GetSnapshot(){return m_Snap;} + + /** Return the gradient direction. idx is 0 based */ + virtual GradientDirectionType GetGradientDirection( unsigned int idx) const + { + if( idx >= m_NumberOfGradientDirections ) + { + itkExceptionMacro( << "Gradient direction " << idx << "does not exist" ); + } + return m_GradientDirectionContainer->ElementAt( idx+1 ); + } + + protected: + DiffusionIntravoxelIncoherentMotionReconstructionImageFilter(); + ~DiffusionIntravoxelIncoherentMotionReconstructionImageFilter() {}; + void PrintSelf(std::ostream& os, Indent indent) const; + + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const + OutputImageRegionType &outputRegionForThread, int); + void AfterThreadedGenerateData(); + + MeasAndBvals ApplyS0Threshold(vnl_vector &meas, vnl_vector &bvals); + + private: + + double myround(double number); + + /** container to hold gradient directions */ + GradientDirectionContainerType::Pointer m_GradientDirectionContainer; + + /** Number of gradient measurements */ + unsigned int m_NumberOfGradientDirections; + + double m_BValue; + + double m_BThres; + + double m_S0Thres; + + double m_DStar; + + IVIM_Method m_Method; + + typename OutputImageType::Pointer m_DMap; + + typename OutputImageType::Pointer m_DStarMap; + +// typename OutputImageType::Pointer m_ADCMap; + + bool m_FitDStar; + + IVIMSnapshot m_Snap; + + bool m_Verbose; + + typename VectorImageType::Pointer m_InternalVectorImage; + + typename InitialFitImageType::Pointer m_InitialFitImage; + + int m_NumberIterations; // for total variation + + vnl_vector m_tmp_allmeas; + + double m_Lambda; + + typename InputImageType::IndexType m_CrossPosition; + +// std::vector m_ADCInds; + + }; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.cpp" +#endif + +#endif //__itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter_h + diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h new file mode 100644 index 0000000000..1e954f24d8 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h @@ -0,0 +1,153 @@ +/*========================================================================= + +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 __itkRegularizedIVIMLocalVariationImageFilter_h +#define __itkRegularizedIVIMLocalVariationImageFilter_h + +#include "itkImageToImageFilter.h" +#include "itkImage.h" + +namespace itk +{ + + template + class IVIMSquaredEuclideanMetric + { + public: + static double Calc(TPixelType p) + { + return p*p; + } + }; + + template<> + class IVIMSquaredEuclideanMetric > + { + public: + static double Calc(itk::Vector p) + { + return p[1]*p[1]; + } + }; + + template<> + class IVIMSquaredEuclideanMetric > + { + public: + static double Calc(itk::VariableLengthVector p) + { + return p.GetSquaredNorm(); + } + }; + + template<> + class IVIMSquaredEuclideanMetric > + { + public: + static double Calc(itk::VariableLengthVector p) + { + return p.GetSquaredNorm(); + } + }; + +/** \class RegularizedIVIMLocalVariationImageFilter + * \brief Calculates the local variation in each pixel + * + * Reference: Tony F. Chan et al., The digital TV filter and nonlinear denoising + * + * \sa Image + * \sa Neighborhood + * \sa NeighborhoodOperator + * \sa NeighborhoodIterator + * + * \ingroup IntensityImageFilters + */ +template +class RegularizedIVIMLocalVariationImageFilter : + public ImageToImageFilter< TInputImage, TOutputImage > +{ +public: + /** Extract dimension from input and output image. */ + itkStaticConstMacro(InputImageDimension, unsigned int, + TInputImage::ImageDimension); + itkStaticConstMacro(OutputImageDimension, unsigned int, + TOutputImage::ImageDimension); + + /** Convenient typedefs for simplifying declarations. */ + typedef TInputImage InputImageType; + typedef TOutputImage OutputImageType; + + /** Standard class typedefs. */ + typedef RegularizedIVIMLocalVariationImageFilter Self; + typedef ImageToImageFilter< InputImageType, OutputImageType> Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(RegularizedIVIMLocalVariationImageFilter, ImageToImageFilter); + + /** Image typedef support. */ + typedef typename InputImageType::PixelType InputPixelType; + typedef typename OutputImageType::PixelType OutputPixelType; + + typedef typename InputImageType::RegionType InputImageRegionType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + + typedef typename InputImageType::SizeType InputSizeType; + + /** MedianImageFilter needs a larger input requested region than + * the output requested region. As such, MedianImageFilter needs + * to provide an implementation for GenerateInputRequestedRegion() + * in order to inform the pipeline execution model. + * + * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ + virtual void GenerateInputRequestedRegion() + throw(InvalidRequestedRegionError); + +protected: + RegularizedIVIMLocalVariationImageFilter(); + virtual ~RegularizedIVIMLocalVariationImageFilter() {} + void PrintSelf(std::ostream& os, Indent indent) const; + + /** MedianImageFilter can be implemented as a multithreaded filter. + * Therefore, this implementation provides a ThreadedGenerateData() + * routine which is called for each processing thread. The output + * image data is allocated automatically by the superclass prior to + * calling ThreadedGenerateData(). ThreadedGenerateData can only + * write to the portion of the output image specified by the + * parameter "outputRegionForThread" + * + * \sa ImageToImageFilter::ThreadedGenerateData(), + * ImageToImageFilter::GenerateData() */ + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int threadId ); + +private: + RegularizedIVIMLocalVariationImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkRegularizedIVIMLocalVariationImageFilter.txx" +#endif + +#endif //RegularizedIVIMLocalVariationImageFilter diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.txx b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.txx new file mode 100644 index 0000000000..a9f2b39aa4 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.txx @@ -0,0 +1,192 @@ +/*========================================================================= + +Program: Insight Segmentation & Registration Toolkit +Language: C++ +Date: $Date: 2006-01-11 19:43:31 $ +Version: $Revision: x $ + +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 +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef _itkRegularizedIVIMLocalVariationImageFilter_txx +#define _itkRegularizedIVIMLocalVariationImageFilter_txx + +#include "itkConstShapedNeighborhoodIterator.h" +#include "itkNeighborhoodInnerProduct.h" +#include "itkImageRegionIterator.h" +#include "itkNeighborhoodAlgorithm.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkOffset.h" +#include "itkProgressReporter.h" +#include "itkVectorImage.h" + +#include +#include + +namespace itk +{ + + template + RegularizedIVIMLocalVariationImageFilter + ::RegularizedIVIMLocalVariationImageFilter() + {} + + template + void + RegularizedIVIMLocalVariationImageFilter + ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError) + { + // call the superclass' implementation of this method + Superclass::GenerateInputRequestedRegion(); + + // get pointers to the input and output + typename Superclass::InputImagePointer inputPtr = + const_cast< TInputImage * >( this->GetInput() ); + typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); + + if ( !inputPtr || !outputPtr ) + { + return; + } + + // get a copy of the input requested region (should equal the output + // requested region) + typename TInputImage::RegionType inputRequestedRegion; + inputRequestedRegion = inputPtr->GetRequestedRegion(); + + // pad the input requested region by 1 + inputRequestedRegion.PadByRadius( 1 ); + + // crop the input requested region at the input's largest possible region + if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) ) + { + inputPtr->SetRequestedRegion( inputRequestedRegion ); + return; + } + else + { + // Couldn't crop the region (requested region is outside the largest + // possible region). Throw an exception. + + // store what we tried to request (prior to trying to crop) + inputPtr->SetRequestedRegion( inputRequestedRegion ); + + // build an exception + InvalidRequestedRegionError e(__FILE__, __LINE__); + e.SetLocation(ITK_LOCATION); + e.SetDescription("Requested region outside possible region."); + e.SetDataObject(inputPtr); + throw e; + } + } + + template< class TInputImage, class TOutputImage> + void + RegularizedIVIMLocalVariationImageFilter< TInputImage, TOutputImage> + ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int threadId) + { + + // Allocate output + typename OutputImageType::Pointer output = this->GetOutput(); + typename InputImageType::ConstPointer input = this->GetInput(); + + itk::Size size; + for( int i=0; i bC; + typename NeighborhoodAlgorithm:: + ImageBoundaryFacesCalculator::FaceListType + faceList = bC(input, outputRegionForThread, size); + + // support progress methods/callbacks + ProgressReporter progress( + this, threadId, outputRegionForThread.GetNumberOfPixels()); + + ZeroFluxNeumannBoundaryCondition nbc; + std::vector pixels; + + // Process each of the boundary faces. These are N-d regions which border + // the edge of the buffer. + for ( typename NeighborhoodAlgorithm:: + ImageBoundaryFacesCalculator::FaceListType::iterator + fit=faceList.begin(); fit != faceList.end(); ++fit) + { + + // iterators over output and input + ImageRegionIterator + output_image_it(output, *fit); + ImageRegionConstIterator + input_image_it(input.GetPointer(), *fit); + + // neighborhood iterator for input image + ConstShapedNeighborhoodIterator + input_image_neighbors_it(size, input, *fit); + typename ConstShapedNeighborhoodIterator:: + OffsetType offset; + input_image_neighbors_it.OverrideBoundaryCondition(&nbc); + input_image_neighbors_it.ClearActiveList(); + for(int i=0; i:: + ConstIterator input_neighbors_it; + for (input_neighbors_it = input_image_neighbors_it.Begin(); + ! input_neighbors_it.IsAtEnd(); + input_neighbors_it++) + { + typename TInputImage::PixelType diffVec = + input_neighbors_it.Get()-input_image_it.Get(); + locVariation += IVIMSquaredEuclideanMetric + ::Calc(diffVec); + } + locVariation = sqrt(locVariation + 0.0001); + output_image_it.Set(locVariation); + + // update iterators + ++input_image_neighbors_it; + ++output_image_it; + ++input_image_it; + + // report progress + progress.CompletedPixel(); + } + } + } + + /** + * Standard "PrintSelf" method + */ + template + void + RegularizedIVIMLocalVariationImageFilter + ::PrintSelf( + std::ostream& os, + Indent indent) const + { + Superclass::PrintSelf( os, indent ); + } + +} // end namespace itk + +#endif //_itkRegularizedIVIMLocalVariationImageFilter_txx diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.h b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.h new file mode 100644 index 0000000000..1d38dd112f --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.h @@ -0,0 +1,132 @@ +/*========================================================================= + +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 __itkRegularizedIVIMReconstructionFilter_h +#define __itkRegularizedIVIMReconstructionFilter_h + +#include "itkRegularizedIVIMReconstructionSingleIteration.h" +#include "itkImageToImageFilter.h" +#include "itkCastImageFilter.h" +#include "itkImage.h" + +namespace itk +{ +/** \class RegularizedIVIMReconstructionFilter + * \brief Applies a total variation denoising filter to an image + * + * Reference: Tony F. Chan et al., The digital TV filter and nonlinear denoising + * + * \sa Image + * \sa Neighborhood + * \sa NeighborhoodOperator + * \sa NeighborhoodIterator + * + * \ingroup IntensityImageFilters + */ +template +class RegularizedIVIMReconstructionFilter : +public ImageToImageFilter< itk::Image, 3>, itk::Image, 3> > +{ +public: + + typedef TInputPixel InputPixelType; + typedef TOutputPixel OutputPixelType; + typedef TRefPixelType RefPixelType; + + /** Convenient typedefs for simplifying declarations. */ + typedef itk::Image, 3> InputImageType; + typedef itk::Image, 3> OutputImageType; + typedef itk::VectorImage RefImageType; + + /** Image typedef support. */ + typedef typename InputImageType::PixelType InputVectorType; + typedef typename OutputImageType::PixelType OutputVectorType; + typedef typename RefImageType::PixelType RefVectorType; + + /** Extract dimension from input and output image. */ + itkStaticConstMacro(InputImageDimension, unsigned int, + InputImageType::ImageDimension); + itkStaticConstMacro(OutputImageDimension, unsigned int, + InputImageType::ImageDimension); + + /** Standard class typedefs. */ + typedef RegularizedIVIMReconstructionFilter Self; + typedef ImageToImageFilter< InputImageType, OutputImageType> Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(RegularizedIVIMReconstructionFilter, ImageToImageFilter); + + typedef typename InputImageType::RegionType InputImageRegionType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + typedef typename RefImageType::RegionType RefImageRegionType; + + typedef typename InputImageType::SizeType InputSizeType; + + typedef RegularizedIVIMReconstructionSingleIteration + + SingleIterationFilterType; + + typedef typename itk::CastImageFilter< InputImageType, OutputImageType > CastType; + + itkSetMacro(Lambda, double); + itkGetMacro(Lambda, double); + + itkSetMacro(NumberIterations, int); + itkGetMacro(NumberIterations, int); + + void SetBValues( vnl_vector bvals ) + { this->m_BValues = bvals; } + vnl_vector GetBValues() + { return this->m_BValues; } + + void SetReferenceImage( typename RefImageType::Pointer refimg ) + { this->m_ReferenceImage = refimg; } + typename RefImageType::Pointer GetReferenceImage() + { return this->m_ReferenceImage; } + +protected: + RegularizedIVIMReconstructionFilter(); + virtual ~RegularizedIVIMReconstructionFilter() {} + void PrintSelf(std::ostream& os, Indent indent) const; + + void GenerateData(); + + double m_Lambda; + + int m_NumberIterations; + + typename RefImageType::Pointer m_ReferenceImage; + + vnl_vector m_BValues; + +private: + RegularizedIVIMReconstructionFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkRegularizedIVIMReconstructionFilter.txx" +#endif + +#endif //__itkRegularizedIVIMReconstructionFilter__ diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.txx b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.txx new file mode 100644 index 0000000000..eedc7ebdb7 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionFilter.txx @@ -0,0 +1,108 @@ +/*========================================================================= + +Program: Insight Segmentation & Registration Toolkit +Language: C++ +Date: $Date: 2006-01-11 19:43:31 $ +Version: $Revision: x $ + +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 +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef _itkRegularizedIVIMReconstructionFilter_txx +#define _itkRegularizedIVIMReconstructionFilter_txx + +#include "itkConstShapedNeighborhoodIterator.h" +#include "itkNeighborhoodInnerProduct.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionConstIterator.h" +#include "itkNeighborhoodAlgorithm.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkOffset.h" +#include "itkProgressReporter.h" + +#include +#include + +namespace itk +{ + + template + RegularizedIVIMReconstructionFilter + ::RegularizedIVIMReconstructionFilter() + { + m_Lambda = 1.0; + } + + + template + void + RegularizedIVIMReconstructionFilter + ::GenerateData() + { + // first we cast the input image to match output type + typename CastType::Pointer infilter = CastType::New(); + infilter->SetInput(this->GetInput()); + infilter->Update(); + typename OutputImageType::Pointer image = infilter->GetOutput(); + + typename SingleIterationFilterType::Pointer filter; + for(int i=0; iSetInput( image.GetPointer() ); + filter->SetOriginalImage( m_ReferenceImage ); + filter->SetBValues(m_BValues); + filter->SetLambda(m_Lambda); + filter->SetNumberOfThreads(this->GetNumberOfThreads()); + filter->UpdateLargestPossibleRegion(); + image = filter->GetOutput(); + std::cout << "Iteration " << i+1 << "/" << + m_NumberIterations << std::endl; + } + + typename OutputImageType::Pointer output = this->GetOutput(); + output->SetOrigin( image->GetOrigin() ); // Set the image origin + output->SetDirection( image->GetDirection() ); // Set the image direction + output->SetSpacing(image->GetSpacing()); + output->SetLargestPossibleRegion( image->GetLargestPossibleRegion() ); + output->SetBufferedRegion( image->GetLargestPossibleRegion() ); + output->Allocate(); + + itk::ImageRegionIterator oit( + output, output->GetLargestPossibleRegion()); + oit.GoToBegin(); + + itk::ImageRegionConstIterator iit( + image, image->GetLargestPossibleRegion()); + iit.GoToBegin(); + + while(!oit.IsAtEnd()) + { + oit.Set(iit.Get()); + ++iit; + ++oit; + } + } + + + /** + * Standard "PrintSelf" method + */ + template + void + RegularizedIVIMReconstructionFilter + ::PrintSelf( + std::ostream& os, + Indent indent) const + { + Superclass::PrintSelf( os, indent ); + } + +} // end namespace itk + +#endif diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h new file mode 100644 index 0000000000..14848f24c9 --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h @@ -0,0 +1,142 @@ +/*========================================================================= + +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 __itkRegularizedIVIMReconstructionSingleIteration_h +#define __itkRegularizedIVIMReconstructionSingleIteration_h + +#include "itkImageToImageFilter.h" +#include "itkImage.h" +#include "itkVectorImage.h" + +namespace itk +{ +/** \class RegularizedIVIMReconstructionSingleIteration + * \brief Applies a total variation denoising filter to an image + * + * Reference: Tony F. Chan et al., The digital TV filter and nonlinear denoising + * + * \sa Image + * \sa Neighborhood + * \sa NeighborhoodOperator + * \sa NeighborhoodIterator + * + * \ingroup IntensityImageFilters + */ +template +class RegularizedIVIMReconstructionSingleIteration : + public ImageToImageFilter< itk::Image, 3>, itk::Image, 3> > +{ +public: + /** Convenient typedefs for simplifying declarations. */ + typedef itk::Image, 3> InputImageType; + typedef itk::Image, 3> OutputImageType; + typedef itk::VectorImage RefImageType; + + /** Extract dimension from input and output image. */ + itkStaticConstMacro(InputImageDimension, unsigned int, + InputImageType::ImageDimension); + itkStaticConstMacro(OutputImageDimension, unsigned int, + OutputImageType::ImageDimension); + + typedef itk::Image LocalVariationImageType; + + /** Standard class typedefs. */ + typedef RegularizedIVIMReconstructionSingleIteration Self; + typedef ImageToImageFilter< InputImageType, OutputImageType> Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(RegularizedIVIMReconstructionSingleIteration, ImageToImageFilter); + + /** Image typedef support. */ + typedef typename InputImageType::PixelType InputPixelType; + typedef typename OutputImageType::PixelType OutputPixelType; + + typedef typename InputImageType::RegionType InputImageRegionType; + typedef typename OutputImageType::RegionType OutputImageRegionType; + + typedef typename InputImageType::SizeType InputSizeType; + + /** A larger input requested region than + * the output requested region is required. + * Therefore, an implementation for GenerateInputRequestedRegion() + * is provided. + * + * \sa ImageToImageFilter::GenerateInputRequestedRegion() */ + virtual void GenerateInputRequestedRegion() + throw(InvalidRequestedRegionError); + + itkSetMacro(Lambda, double); + itkGetMacro(Lambda, double); + + void SetBValues(vnl_vector bvals) + { this->m_BValues = bvals; } + vnl_vector GetBValues() + { return this->m_BValues; } + + void SetOriginalImage(RefImageType* in) + { this->m_OriginalImage = in; } + typename RefImageType::Pointer GetOriginialImage() + { return this->m_OriginalImage; } + +protected: + RegularizedIVIMReconstructionSingleIteration(); + virtual ~RegularizedIVIMReconstructionSingleIteration() {} + void PrintSelf(std::ostream& os, Indent indent) const; + + /** MedianImageFilter can be implemented as a multithreaded filter. + * Therefore, this implementation provides a ThreadedGenerateData() + * routine which is called for each processing thread. The output + * image data is allocated automatically by the superclass prior to + * calling ThreadedGenerateData(). ThreadedGenerateData can only + * write to the portion of the output image specified by the + * parameter "outputRegionForThread" + * + * \sa ImageToImageFilter::ThreadedGenerateData(), + * ImageToImageFilter::GenerateData() */ + void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int threadId ); + + void BeforeThreadedGenerateData(); + + typename LocalVariationImageType::Pointer m_LocalVariation; + + typename RefImageType::Pointer m_OriginalImage; + + double m_Lambda; + + vnl_vector m_BValues; + +private: + + RegularizedIVIMReconstructionSingleIteration(const Self&); + + void operator=(const Self&); + +}; + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkRegularizedIVIMReconstructionSingleIteration.txx" +#endif + +#endif //__itkRegularizedIVIMReconstructionSingleIteration__ diff --git a/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.txx b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.txx new file mode 100644 index 0000000000..5564ab77fd --- /dev/null +++ b/Modules/DiffusionImaging/Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.txx @@ -0,0 +1,310 @@ +/*========================================================================= + +Program: Insight Segmentation & Registration Toolkit +Language: C++ +Date: $Date: 2006-01-11 19:43:31 $ +Version: $Revision: x $ + +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 +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#ifndef _itkRegularizedIVIMReconstructionSingleIteration_txx +#define _itkRegularizedIVIMReconstructionSingleIteration_txx + +// itk includes +#include "itkConstShapedNeighborhoodIterator.h" +#include "itkNeighborhoodInnerProduct.h" +#include "itkImageRegionIterator.h" +#include "itkNeighborhoodAlgorithm.h" +#include "itkZeroFluxNeumannBoundaryCondition.h" +#include "itkOffset.h" +#include "itkProgressReporter.h" +#include "itkRegularizedIVIMLocalVariationImageFilter.h" + +// other includes +#include +#include + +#define IVIM_FOO -100000 + +namespace itk +{ + + /** + * constructor + */ + template + RegularizedIVIMReconstructionSingleIteration + ::RegularizedIVIMReconstructionSingleIteration() + { + m_Lambda = 1.0; + m_LocalVariation = LocalVariationImageType::New(); + } + + /** + * generate requested region + */ + template + void + RegularizedIVIMReconstructionSingleIteration + ::GenerateInputRequestedRegion() throw (InvalidRequestedRegionError) + { + // call the superclass' implementation of this method + Superclass::GenerateInputRequestedRegion(); + + // get pointers to the input and output + typename Superclass::InputImagePointer inputPtr = + const_cast< InputImageType * >( this->GetInput() ); + typename Superclass::OutputImagePointer outputPtr = this->GetOutput(); + + if ( !inputPtr || !outputPtr ) + { + return; + } + + // get a copy of the input requested region (should equal the output + // requested region) + typename InputImageType::RegionType inputRequestedRegion; + inputRequestedRegion = inputPtr->GetRequestedRegion(); + + // pad the input requested region by 1 + inputRequestedRegion.PadByRadius( 1 ); + + // crop the input requested region at the input's largest possible region + if ( inputRequestedRegion.Crop(inputPtr->GetLargestPossibleRegion()) ) + { + inputPtr->SetRequestedRegion( inputRequestedRegion ); + return; + } + else + { + // Couldn't crop the region (requested region is outside the largest + // possible region). Throw an exception. + + // store what we tried to request (prior to trying to crop) + inputPtr->SetRequestedRegion( inputRequestedRegion ); + + // build an exception + InvalidRequestedRegionError e(__FILE__, __LINE__); + e.SetLocation(ITK_LOCATION); + e.SetDescription("Requested region outside possible region."); + e.SetDataObject(inputPtr); + throw e; + } + } + + + /** + * generate output + */ + template + void + RegularizedIVIMReconstructionSingleIteration + ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int threadId) + { + + typename OutputImageType::Pointer output = this->GetOutput(); + typename InputImageType::ConstPointer input = this->GetInput(); + + // Find the data-set boundary "faces" + itk::Size size; + for( int i=0; i bC; + typename NeighborhoodAlgorithm:: + ImageBoundaryFacesCalculator::FaceListType + faceList = bC(input, outputRegionForThread, size); + + NeighborhoodAlgorithm:: + ImageBoundaryFacesCalculator lv_bC; + typename NeighborhoodAlgorithm:: + ImageBoundaryFacesCalculator::FaceListType + lv_faceList = lv_bC(m_LocalVariation, outputRegionForThread, size); + + ZeroFluxNeumannBoundaryCondition nbc; + ZeroFluxNeumannBoundaryCondition lv_nbc; + std::vector ws; + std::vector hs; + + typename NeighborhoodAlgorithm:: + ImageBoundaryFacesCalculator::FaceListType::iterator + lv_fit=lv_faceList.begin(); + + // Process each of the boundary faces. These are N-d regions which border + // the edge of the buffer. + for ( typename NeighborhoodAlgorithm:: + ImageBoundaryFacesCalculator::FaceListType::iterator + fit=faceList.begin(); fit != faceList.end(); ++fit) + { + + // iterators over output, input, original and local variation image + ImageRegionIterator output_image_it = + ImageRegionIterator(output, *fit); + ImageRegionConstIterator input_image_it = + ImageRegionConstIterator(input, *fit); + ImageRegionConstIterator orig_image_it = + ImageRegionConstIterator(m_OriginalImage, *fit); + ImageRegionConstIterator loc_var_image_it = + ImageRegionConstIterator( + m_LocalVariation, *fit); + + // neighborhood in input image + ConstShapedNeighborhoodIterator + input_image_neighbors_it(size, input, *fit); + typename ConstShapedNeighborhoodIterator:: + OffsetType offset; + input_image_neighbors_it.OverrideBoundaryCondition(&nbc); + input_image_neighbors_it.ClearActiveList(); + for(int i=0; i + loc_var_image_neighbors_it(size, m_LocalVariation, *lv_fit); + loc_var_image_neighbors_it.OverrideBoundaryCondition(&lv_nbc); + loc_var_image_neighbors_it.ClearActiveList(); + for(int i=0; i:: + ConstIterator loc_var_neighbors_it; + for (loc_var_neighbors_it = loc_var_image_neighbors_it.Begin(); + ! loc_var_neighbors_it.IsAtEnd(); + loc_var_neighbors_it++) + { + // w_alphabeta(u) = + // 1 / ||nabla_alpha(u)||_a + 1 / ||nabla_beta(u)||_a + ws[count] = + locvar_alpha_inv + (1.0/(double)loc_var_neighbors_it.Get()); + wsum += ws[count++]; + + //MITK_INFO << "nb var: " << count << ": " << loc_var_neighbors_it.Get(); + } + //MITK_INFO << "wsum: " << wsum; + + // h_alphaalpha * u_alpha^zero + typename RefImageType::PixelType orig = orig_image_it.Get(); +// vnl_vector estim(orig.GetSize()); +// vnl_matrix diff(orig.GetSize(),1); +// vnl_matrix estimdash(orig.GetSize(),2); + vnl_vector_fixed step; + step[0] = 0; step[1]=0; + for(int ind=0; ind:: + ConstIterator input_neighbors_it; + for (input_neighbors_it = input_image_neighbors_it.Begin(); + ! input_neighbors_it.IsAtEnd(); + input_neighbors_it++) + { + step[1] += (input_neighbors_it.Get()[1] - input_image_it.Get()[1]) * (ws[count++] / (m_Lambda+wsum)); + } + + //MITK_INFO << "stepfinal: " << step[0] << "; " << step[1]; + + // set output result + OutputPixelType out; + out[0] = input_image_it.Get()[0] + .001*step[0]; + out[1] = input_image_it.Get()[1] + .00001*step[1]; + output_image_it.Set( out ); + + //MITK_INFO << "(" << input_image_it.Get()[0] << " ; " << input_image_it.Get()[1] << ") => (" << out[0] << " ; " << out[1] << ")"; + // increment iterators + ++output_image_it; + ++input_image_it; + ++orig_image_it; + ++loc_var_image_it; + ++input_image_neighbors_it; + ++loc_var_image_neighbors_it; + + } + + ++lv_fit; + } + } + + /** + * first calculate local variation in the image + */ + template + void + RegularizedIVIMReconstructionSingleIteration + ::BeforeThreadedGenerateData() + { + typedef typename itk::RegularizedIVIMLocalVariationImageFilter + FilterType; + typename FilterType::Pointer filter = FilterType::New(); + filter->SetInput(this->GetInput(0)); + filter->SetNumberOfThreads(this->GetNumberOfThreads()); + filter->Update(); + this->m_LocalVariation = filter->GetOutput(); + } + + /** + * Standard "PrintSelf" method + */ + template + void + RegularizedIVIMReconstructionSingleIteration + ::PrintSelf( + std::ostream& os, + Indent indent) const + { + Superclass::PrintSelf( os, indent ); + } + +} // end namespace itk + +#endif diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 1e86e59231..9a079ea087 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,156 +1,160 @@ SET(CPP_FILES # 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 IODataStructures/FiberBundle/mitkParticleGrid.cpp # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures -> Tbss IODataStructures/TbssImages/mitkTbssImageSource.cpp IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp IODataStructures/TbssImages/mitkTbssImage.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleMapper2D.cpp Rendering/mitkFiberBundleMapper3D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp ) SET(H_FILES # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleMapper2D.h Rendering/mitkFiberBundleMapper3D.h Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h + Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h + Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h + Reconstruction/itkRegularizedIVIMReconstructionFilter.h + Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/FiberBundle/itkSlowPolyLineParametricPath.h IODataStructures/TbssImages/mitkTbssImage.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h # Tractography Tractography/itkGibbsTrackingFilter.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 Algorithms/itkTractsToDWIImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h ) SET( TOOL_FILES ) IF(WIN32) ENDIF(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp )