diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake index 42a2b36f57..90a9f5c50a 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake @@ -1,86 +1,91 @@ SET(SRC_CPP_FILES QmitkODFDetailsWidget.cpp QmitkODFRenderWidget.cpp + QmitkPartialVolumeAnalysisWidget.cpp ) SET(INTERNAL_CPP_FILES mitkPluginActivator.cpp QmitkQBallReconstructionView.cpp QmitkPreprocessingView.cpp QmitkDiffusionDicomImportView.cpp QmitkDiffusionQuantificationView.cpp QmitkTensorReconstructionView.cpp QmitkDiffusionImagingPublicPerspective.cpp QmitkControlVisualizationPropertiesView.cpp QmitkODFDetailsView.cpp QmitkGlobalFiberTrackingView.cpp QmitkFiberBundleOperationsView.cpp QmitkFiberBundleDeveloperView.cpp + QmitkPartialVolumeAnalysisView.cpp ) SET(UI_FILES src/internal/QmitkQBallReconstructionViewControls.ui src/internal/QmitkPreprocessingViewControls.ui src/internal/QmitkDiffusionDicomImportViewControls.ui src/internal/QmitkDiffusionQuantificationViewControls.ui src/internal/QmitkTensorReconstructionViewControls.ui src/internal/QmitkControlVisualizationPropertiesViewControls.ui src/internal/QmitkODFDetailsViewControls.ui src/internal/QmitkGlobalFiberTrackingViewControls.ui src/internal/QmitkFiberBundleOperationsViewControls.ui src/internal/QmitkFiberBundleDeveloperViewControls.ui + src/internal/QmitkPartialVolumeAnalysisViewControls.ui ) SET(MOC_H_FILES src/internal/mitkPluginActivator.h src/internal/QmitkQBallReconstructionView.h src/internal/QmitkPreprocessingView.h src/internal/QmitkDiffusionDicomImportView.h src/internal/QmitkDiffusionImagingPublicPerspective.h src/internal/QmitkDiffusionQuantificationView.h src/internal/QmitkTensorReconstructionView.h src/internal/QmitkControlVisualizationPropertiesView.h src/internal/QmitkODFDetailsView.h src/QmitkODFRenderWidget.h src/QmitkODFDetailsWidget.h src/internal/QmitkGlobalFiberTrackingView.h src/internal/QmitkFiberBundleOperationsView.h src/internal/QmitkFiberBundleDeveloperView.h + src/internal/QmitkPartialVolumeAnalysisView.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/GlobalTracking.png resources/FiberBundleOperations.png + resources/PartialVolumeAnalysis_24.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 1607fa9d65..0f57022457 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml @@ -1,84 +1,89 @@ - + - + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_24.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_24.png new file mode 100644 index 0000000000..1fee2ecf3e Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_24.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_48.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_48.png new file mode 100644 index 0000000000..a664faeb2c Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_48.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_48.xcf b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_48.xcf new file mode 100644 index 0000000000..24328a1355 Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_48.xcf differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_64.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_64.png new file mode 100644 index 0000000000..bf1e4b94a5 Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/PartialVolumeAnalysis_64.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkPartialVolumeAnalysisWidget.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkPartialVolumeAnalysisWidget.cpp new file mode 100644 index 0000000000..f362604cd0 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkPartialVolumeAnalysisWidget.cpp @@ -0,0 +1,127 @@ +/*========================================================================= +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 "QmitkPartialVolumeAnalysisWidget.h" + +#include "mitkHistogramGenerator.h" +#include "mitkPartialVolumeAnalysisClusteringCalculator.h" + +#include +#include +#include + + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +//#include + +QmitkPartialVolumeAnalysisWidget::QmitkPartialVolumeAnalysisWidget( QWidget * parent ) + : QmitkPlotWidget(parent) +{ +// this->SetAxisTitle( QwtPlot::xBottom, "Grayvalue" ); + // this->SetAxisTitle( QwtPlot::yLeft, "Probability" ); +// this->Replot(); + m_Plot->setCanvasLineWidth(0); + m_Plot->setMargin(0); +} + + + +QmitkPartialVolumeAnalysisWidget::~QmitkPartialVolumeAnalysisWidget() +{ +} + + +void QmitkPartialVolumeAnalysisWidget::DrawGauss() +{ + +} + + +void QmitkPartialVolumeAnalysisWidget::ClearItemModel() +{ + +} + +void QmitkPartialVolumeAnalysisWidget::SetParameters( ParamsType *params, ResultsType *results, HistType *hist ) +{ + this->Clear(); + + if(params != 0 && results != 0) + { + params->Print(); + + for(unsigned int i=0; iGetXVals()); + m_Vals.push_back(hist->GetHVals()); + + QPen pen( Qt::SolidLine ); + pen.setWidth(2); + + pen.setColor(Qt::black); + int curveId = this->InsertCurve( "histogram" ); + this->SetCurveData( curveId, (*m_Vals[0]), (*m_Vals[1]) ); + this->SetCurvePen( curveId, pen ); + // this->SetCurveTitle( curveId, "Image Histogram" ); + + std::vector *fiberVals = new std::vector(results->GetFiberVals()); + curveId = this->InsertCurve( "fiber" ); + this->SetCurveData( curveId, (*hist->GetXVals()), (*fiberVals) ); + this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); + this->SetCurveBrush(curveId, QBrush(QColor::fromRgbF(0,1,0,.5), Qt::SolidPattern)); + m_Vals.push_back(fiberVals); + + std::vector *nonFiberVals = new std::vector(results->GetNonFiberVals()); + curveId = this->InsertCurve( "nonfiber" ); + this->SetCurveData( curveId, (*hist->GetXVals()), (*nonFiberVals) ); + this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); + this->SetCurveBrush(curveId, QBrush(QColor::fromRgbF(1,0,0,.5), Qt::SolidPattern)); + m_Vals.push_back(nonFiberVals); + + std::vector *mixedVals = new std::vector(results->GetMixedVals()); + curveId = this->InsertCurve( "mixed" ); + this->SetCurveData( curveId, (*hist->GetXVals()), (*mixedVals) ); + this->SetCurvePen( curveId, QPen( Qt::NoPen ) ); + this->SetCurveBrush(curveId, QBrush(QColor::fromRgbF(.7,.7,.7,.5), Qt::SolidPattern)); + m_Vals.push_back(mixedVals); + + pen.setColor(Qt::blue); + std::vector *combiVals = new std::vector(results->GetCombiVals()); + curveId = this->InsertCurve( "combi" ); + this->SetCurveData( curveId, (*hist->GetXVals()), (*combiVals) ); + this->SetCurvePen( curveId, pen ); + m_Vals.push_back(combiVals); + + + + } + + this->Replot(); + +} + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkPartialVolumeAnalysisWidget.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkPartialVolumeAnalysisWidget.h new file mode 100644 index 0000000000..dc6711e9d0 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkPartialVolumeAnalysisWidget.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 QmitkPartialVolumeAnalysisWidget_H_ +#define QmitkPartialVolumeAnalysisWidget_H_ + +#include "QmitkPlotWidget.h" +#include + +#include "QmitkHistogram.h" +#include "QmitkExtExports.h" +#include "mitkImage.h" +#include "mitkPlanarFigure.h" +#include "mitkPartialVolumeAnalysisClusteringCalculator.h" + +#include +#include +#include + +#include + +#include +#include + +#include + + + + +/** + * \brief Widget for displaying image histograms based on the vtkQtChart + * framework + */ +class MitkDiffusionImaging_EXPORT QmitkPartialVolumeAnalysisWidget : public QmitkPlotWidget +{ + +public: + typedef mitk::Image::HistogramType HistogramType; + typedef mitk::Image::HistogramType::ConstIterator HistogramConstIteratorType; + + typedef mitk::PartialVolumeAnalysisClusteringCalculator ClusteringType; + typedef ClusteringType::ParamsType ParamsType; + typedef ClusteringType::ClusterResultType ResultsType; + typedef ClusteringType::HistType HistType; + + /** \brief Set histogram to be displayed directly. */ + void SetParameters( ParamsType *params, ResultsType *results, HistType *hist ); + + QmitkPartialVolumeAnalysisWidget( QWidget * /*parent = 0 */); + virtual ~QmitkPartialVolumeAnalysisWidget(); + + void DrawGauss(); + + void ClearItemModel(); + + std::vector< std::vector* > m_Vals; +}; + +#endif /* QmitkPartialVolumeAnalysisWidget_H_ */ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisView.cpp new file mode 100644 index 0000000000..08ed7f21f7 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisView.cpp @@ -0,0 +1,2044 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-05-22 11:00:35 +0200 (Fr, 22 Mai 2009) $ +Version: $Revision: 10185 $ + +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 "QmitkPartialVolumeAnalysisView.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +#include "QmitkStdMultiWidget.h" +#include "QmitkSliderNavigatorWidget.h" + +#include "mitkNodePredicateDataType.h" +#include "mitkNodePredicateOr.h" +#include "mitkImageTimeSelector.h" +#include "mitkProperties.h" +#include "mitkProgressBar.h" + +// Includes for image processing +#include "mitkImageCast.h" +#include "mitkImageToItk.h" +#include "mitkITKImageImport.h" +#include "mitkDataNodeObject.h" +#include "mitkNodePredicateData.h" + +#include "mitkPlanarFigureInteractor.h" +#include "mitkGlobalInteraction.h" +#include "mitkTensorImage.h" + +#include "mitkPlanarCircle.h" +#include "mitkPlanarRectangle.h" +#include "mitkPlanarPolygon.h" + +#include "mitkPartialVolumeAnalysisClusteringCalculator.h" + +#include +#include "itkTensorDerivedMeasurementsFilter.h" +#include "itkDiffusionTensor3D.h" +#include "itkCartesianToPolarVectorImageFilter.h" +#include "itkPolarToCartesianVectorImageFilter.h" +#include "itkBinaryThresholdImageFilter.h" +#include "itkMaskImageFilter.h" + +#include "itkImageMomentsCalculator.h" + +#include + +#define PVA_PI 3.141592653589793238462643383 + +class QmitkRequestStatisticsUpdateEvent : public QEvent +{ +public: + enum Type + { + StatisticsUpdateRequest = QEvent::MaxUser - 1025 + }; + + QmitkRequestStatisticsUpdateEvent() + : QEvent( (QEvent::Type) StatisticsUpdateRequest ) {}; +}; + + + +typedef itk::Image ImageType; +typedef itk::Image FloatImageType; +typedef itk::Image, 3> VectorImageType; + +inline bool my_isnan(float x) +{ + volatile float d = x; + + if(d!=d) + return true; + + if(d==d) + return false; + return d != d; + +} + +QmitkPartialVolumeAnalysisView::QmitkPartialVolumeAnalysisView(QObject */*parent*/, const char */*name*/) + : QmitkFunctionality(), + m_Controls( NULL ), + m_TimeStepperAdapter( NULL ), + m_MeasurementInfoRenderer(0), + m_MeasurementInfoAnnotation(0), + m_SelectedImageNodes( ), + m_SelectedImage( NULL ), + m_SelectedMaskNode( NULL ), + m_SelectedImageMask( NULL ), + m_SelectedPlanarFigureNodes(0), + m_SelectedPlanarFigure( NULL ), + m_IsTensorImage(false), + m_FAImage(0), + m_RDImage(0), + m_ADImage(0), + m_MDImage(0), + m_CAImage(0), + // m_DirectionImage(0), + m_DirectionComp1Image(0), + m_DirectionComp2Image(0), + m_AngularErrorImage(0), + m_SelectedRenderWindow(NULL), + m_LastRenderWindow(NULL), + m_ImageObserverTag( -1 ), + m_ImageMaskObserverTag( -1 ), + m_PlanarFigureObserverTag( -1 ), + m_CurrentStatisticsValid( false ), + m_StatisticsUpdatePending( false ), + m_GaussianSigmaChangedSliding(false), + m_NumberBinsSliding(false), + m_UpsamplingChangedSliding(false), + m_ClusteringResult(NULL), + m_EllipseCounter(0), + m_RectangleCounter(0), + m_PolygonCounter(0), + m_CurrentFigureNodeInitialized(false), + m_QuantifyClass(2), + m_IconTexOFF(new QIcon(":/QmitkPartialVolumeAnalysisView/texIntOFFIcon.png")), + m_IconTexON(new QIcon(":/QmitkPartialVolumeAnalysisView/texIntONIcon.png")), + m_TexIsOn(true) +{ + +} + + +QmitkPartialVolumeAnalysisView::~QmitkPartialVolumeAnalysisView() +{ + if ( m_SelectedImage.IsNotNull() ) + m_SelectedImage->RemoveObserver( m_ImageObserverTag ); + if ( m_SelectedImageMask.IsNotNull() ) + m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); + if ( m_SelectedPlanarFigure.IsNotNull() ) + { + m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); + m_SelectedPlanarFigure->RemoveObserver( m_InitializedObserverTag ); + } + + this->GetDefaultDataStorage()->AddNodeEvent -= mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeAddedInDataStorage ); + + m_SelectedPlanarFigureNodes->NodeChanged.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); + + m_SelectedPlanarFigureNodes->NodeRemoved.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); + + m_SelectedPlanarFigureNodes->PropertyChanged.RemoveListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); + + m_SelectedImageNodes->NodeChanged.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); + + m_SelectedImageNodes->NodeRemoved.RemoveListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); + + m_SelectedImageNodes->PropertyChanged.RemoveListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); +} + + +void QmitkPartialVolumeAnalysisView::CreateQtPartControl(QWidget *parent) +{ + if (m_Controls == NULL) + { + m_Controls = new Ui::QmitkPartialVolumeAnalysisViewControls; + m_Controls->setupUi(parent); + this->CreateConnections(); + + m_Controls->m_ErrorMessageLabel->hide(); + } + + SetHistogramVisibility(); + + m_Controls->m_TextureIntON->setIcon(*m_IconTexON); + + m_Controls->m_SimilarAnglesFrame->setVisible(false); + m_Controls->m_SimilarAnglesLabel->setVisible(false); + + vtkTextProperty *textProp = vtkTextProperty::New(); + textProp->SetColor(1.0, 1.0, 1.0); + + m_MeasurementInfoAnnotation = vtkCornerAnnotation::New(); + m_MeasurementInfoAnnotation->SetMaximumFontSize(12); + m_MeasurementInfoAnnotation->SetTextProperty(textProp); + + m_MeasurementInfoRenderer = vtkRenderer::New(); + m_MeasurementInfoRenderer->AddActor(m_MeasurementInfoAnnotation); + + m_SelectedPlanarFigureNodes = mitk::DataStorageSelection::New(this->GetDefaultDataStorage(), false); + + m_SelectedPlanarFigureNodes->NodeChanged.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); + + m_SelectedPlanarFigureNodes->NodeRemoved.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); + + m_SelectedPlanarFigureNodes->PropertyChanged.AddListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); + + m_SelectedImageNodes = mitk::DataStorageSelection::New(this->GetDefaultDataStorage(), false); + + m_SelectedImageNodes->PropertyChanged.AddListener( mitk::MessageDelegate2( this, &QmitkPartialVolumeAnalysisView::PropertyChanged ) ); + + m_SelectedImageNodes->NodeChanged.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeChanged ) ); + + m_SelectedImageNodes->NodeRemoved.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeRemoved ) ); + + this->GetDefaultDataStorage()->AddNodeEvent.AddListener( mitk::MessageDelegate1( this, &QmitkPartialVolumeAnalysisView::NodeAddedInDataStorage ) ); + + Select(NULL,true,true); + + SetAdvancedVisibility(); +} + +void QmitkPartialVolumeAnalysisView::SetHistogramVisibility() +{ + m_Controls->m_HistogramWidget->setVisible(m_Controls->m_DisplayHistogramCheckbox->isChecked()); +} + +void QmitkPartialVolumeAnalysisView::SetAdvancedVisibility() +{ + m_Controls->frame_7->setVisible(m_Controls->m_AdvancedCheckbox->isChecked()); +} + +void QmitkPartialVolumeAnalysisView::CreateConnections() +{ + if ( m_Controls ) + { + + connect( m_Controls->m_DisplayHistogramCheckbox, SIGNAL( clicked() ) + , this, SLOT( SetHistogramVisibility() ) ); + + connect( m_Controls->m_AdvancedCheckbox, SIGNAL( clicked() ) + , this, SLOT( SetAdvancedVisibility() ) ); + + + connect( m_Controls->m_NumberBinsSlider, SIGNAL( sliderReleased () ), + this, SLOT( NumberBinsReleasedSlider( ) ) ); + connect( m_Controls->m_UpsamplingSlider, SIGNAL( sliderReleased( ) ), + this, SLOT( UpsamplingReleasedSlider( ) ) ); + connect( m_Controls->m_GaussianSigmaSlider, SIGNAL( sliderReleased( ) ), + this, SLOT( GaussianSigmaReleasedSlider( ) ) ); + connect( m_Controls->m_SimilarAnglesSlider, SIGNAL( sliderReleased( ) ), + this, SLOT( SimilarAnglesReleasedSlider( ) ) ); + + connect( m_Controls->m_NumberBinsSlider, SIGNAL( valueChanged (int) ), + this, SLOT( NumberBinsChangedSlider( int ) ) ); + connect( m_Controls->m_UpsamplingSlider, SIGNAL( valueChanged( int ) ), + this, SLOT( UpsamplingChangedSlider( int ) ) ); + connect( m_Controls->m_GaussianSigmaSlider, SIGNAL( valueChanged( int ) ), + this, SLOT( GaussianSigmaChangedSlider( int ) ) ); + connect( m_Controls->m_SimilarAnglesSlider, SIGNAL( valueChanged( int ) ), + this, SLOT( SimilarAnglesChangedSlider(int) ) ); + + connect( m_Controls->m_OpacitySlider, SIGNAL( valueChanged( int ) ), + this, SLOT( OpacityChangedSlider(int) ) ); + + connect( (QObject*)(m_Controls->m_ButtonCopyHistogramToClipboard), SIGNAL(clicked()),(QObject*) this, SLOT(ToClipBoard())); + + connect( m_Controls->m_CircleButton, SIGNAL( clicked() ) + , this, SLOT( ActionDrawEllipseTriggered() ) ); + + connect( m_Controls->m_RectangleButton, SIGNAL( clicked() ) + , this, SLOT( ActionDrawRectangleTriggered() ) ); + + connect( m_Controls->m_PolygonButton, SIGNAL( clicked() ) + , this, SLOT( ActionDrawPolygonTriggered() ) ); + + connect( m_Controls->m_GreenRadio, SIGNAL( clicked(bool) ) + , this, SLOT( GreenRadio(bool) ) ); + + connect( m_Controls->m_PartialVolumeRadio, SIGNAL( clicked(bool) ) + , this, SLOT( PartialVolumeRadio(bool) ) ); + + connect( m_Controls->m_BlueRadio, SIGNAL( clicked(bool) ) + , this, SLOT( BlueRadio(bool) ) ); + + connect( m_Controls->m_AllRadio, SIGNAL( clicked(bool) ) + , this, SLOT( AllRadio(bool) ) ); + + connect( m_Controls->m_EstimateCircle, SIGNAL( clicked() ) + , this, SLOT( EstimateCircle() ) ); + + connect( (QObject*)(m_Controls->m_TextureIntON), SIGNAL(clicked()), this, SLOT(TextIntON()) ); + + } +} + +void QmitkPartialVolumeAnalysisView::EstimateCircle() +{ + typedef itk::Image SegImageType; + SegImageType::Pointer mask_itk = SegImageType::New(); + + typedef mitk::ImageToItk CastType; + CastType::Pointer caster = CastType::New(); + caster->SetInput(m_SelectedImageMask); + caster->Update(); + + typedef itk::ImageMomentsCalculator< SegImageType > MomentsType; + MomentsType::Pointer momentsCalc = MomentsType::New(); + momentsCalc->SetImage(caster->GetOutput()); + momentsCalc->Compute(); + MomentsType::VectorType cog = momentsCalc->GetCenterOfGravity(); + MomentsType::MatrixType axes = momentsCalc->GetPrincipalAxes(); + MomentsType::VectorType moments = momentsCalc->GetPrincipalMoments(); + + // moments-coord conversion + // third coordinate min oder max? + + // max-min = extent + + MomentsType::AffineTransformPointer trafo = momentsCalc->GetPhysicalAxesToPrincipalAxesTransform(); + + itk::ImageRegionIterator + itimage(caster->GetOutput(), caster->GetOutput()->GetLargestPossibleRegion()); + itimage = itimage.Begin(); + + double max = -9999999999.0; + double min = 9999999999.0; + + while( !itimage.IsAtEnd() ) + { + if(itimage.Get()) + { + ImageType::IndexType index = itimage.GetIndex(); + + itk::Point point; + caster->GetOutput()->TransformIndexToPhysicalPoint(index,point); + + itk::Point newPoint; + newPoint = trafo->TransformPoint(point); + + if(newPoint[2]max) + max = newPoint[2]; + } + ++itimage; + } + + double extent = max - min; + MITK_INFO << "EXTENT = " << extent; + + mitk::Point3D origin; + mitk::Vector3D right, bottom, normal; + + double factor = 1000.0; + mitk::FillVector3D(origin, cog[0]-factor*axes[1][0]-factor*axes[2][0], + cog[1]-factor*axes[1][1]-factor*axes[2][1], + cog[2]-factor*axes[1][2]-factor*axes[2][2]); + // mitk::FillVector3D(normal, axis[0][0],axis[0][1],axis[0][2]); + mitk::FillVector3D(bottom, 2*factor*axes[1][0], 2*factor*axes[1][1], 2*factor*axes[1][2]); + mitk::FillVector3D(right, 2*factor*axes[2][0], 2*factor*axes[2][1], 2*factor*axes[2][2]); + + mitk::PlaneGeometry::Pointer planegeometry = mitk::PlaneGeometry::New(); + planegeometry->InitializeStandardPlane(right.Get_vnl_vector(), bottom.Get_vnl_vector()); + planegeometry->SetOrigin(origin); + + double len1 = sqrt(axes[1][0]*axes[1][0] + axes[1][1]*axes[1][1] + axes[1][2]*axes[1][2]); + double len2 = sqrt(axes[2][0]*axes[2][0] + axes[2][1]*axes[2][1] + axes[2][2]*axes[2][2]); + + mitk::Point2D point1; + point1[0] = factor*len1; + point1[1] = factor*len2; + + mitk::Point2D point2; + point2[0] = factor*len1+extent*.5; + point2[1] = factor*len2; + + mitk::PlanarCircle::Pointer circle = mitk::PlanarCircle::New(); + circle->SetGeometry2D(planegeometry); + circle->PlaceFigure( point1 ); + circle->SetControlPoint(0,point1); + circle->SetControlPoint(1,point2); + //circle->SetCurrentControlPoint( point2 ); + + mitk::PlanarFigure::PolyLineType polyline = circle->GetPolyLine( 0 ); + MITK_INFO << "SIZE of planar figure polyline: " << polyline.size(); + + AddFigureToDataStorage(circle, "Circle"); + +} + +void QmitkPartialVolumeAnalysisView::StdMultiWidgetAvailable( QmitkStdMultiWidget& stdMultiWidget ) +{ + QmitkFunctionality::StdMultiWidgetAvailable(stdMultiWidget); +} + +bool QmitkPartialVolumeAnalysisView::AssertDrawingIsPossible(bool checked) +{ + if (m_SelectedImageNodes->GetNode().IsNull()) + { + checked = false; + this->HandleException("Please select an image!", this->m_Parent, true); + return false; + } + + //this->GetActiveStdMultiWidget()->SetWidgetPlanesVisibility(false); + + return checked; +} + +void QmitkPartialVolumeAnalysisView::ActionDrawEllipseTriggered() +{ + bool checked = m_Controls->m_CircleButton->isChecked(); + if(!this->AssertDrawingIsPossible(checked)) + return; + + mitk::PlanarCircle::Pointer figure = mitk::PlanarCircle::New(); + this->AddFigureToDataStorage(figure, QString("Circle%1").arg(++m_EllipseCounter)); + + MITK_INFO << "PlanarCircle created ..."; +} + +void QmitkPartialVolumeAnalysisView::ActionDrawRectangleTriggered() +{ + bool checked = m_Controls->m_RectangleButton->isChecked(); + if(!this->AssertDrawingIsPossible(checked)) + return; + + mitk::PlanarRectangle::Pointer figure = mitk::PlanarRectangle::New(); + this->AddFigureToDataStorage(figure, QString("Rectangle%1").arg(++m_RectangleCounter)); + + MITK_INFO << "PlanarRectangle created ..."; +} + +void QmitkPartialVolumeAnalysisView::ActionDrawPolygonTriggered() +{ + bool checked = m_Controls->m_PolygonButton->isChecked(); + if(!this->AssertDrawingIsPossible(checked)) + return; + + mitk::PlanarPolygon::Pointer figure = mitk::PlanarPolygon::New(); + figure->ClosedOn(); + this->AddFigureToDataStorage(figure, QString("Polygon%1").arg(++m_PolygonCounter)); + + MITK_INFO << "PlanarPolygon created ..."; +} + +void QmitkPartialVolumeAnalysisView::AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, + const char *propertyKey, mitk::BaseProperty *property ) +{ + + mitk::DataNode::Pointer newNode = mitk::DataNode::New(); + newNode->SetName(name.toStdString()); + newNode->SetData(figure); + + // Add custom property, if available + if ( (propertyKey != NULL) && (property != NULL) ) + { + newNode->AddProperty( propertyKey, property ); + } + + // figure drawn on the topmost layer / image + this->GetDataStorage()->Add(newNode, m_SelectedImageNodes->GetNode() ); + + std::vector selectedNodes = GetDataManagerSelection(); + for(unsigned int i = 0; i < selectedNodes.size(); i++) + { + selectedNodes[i]->SetSelected(false); + } + selectedNodes = m_SelectedPlanarFigureNodes->GetNodes(); + for(unsigned int i = 0; i < selectedNodes.size(); i++) + { + selectedNodes[i]->SetSelected(false); + } + newNode->SetSelected(true); + + Select(newNode); +} + +void QmitkPartialVolumeAnalysisView::PlanarFigureInitialized() +{ + if(m_SelectedPlanarFigureNodes->GetNode().IsNull()) + return; + + m_CurrentFigureNodeInitialized = true; + + this->Select(m_SelectedPlanarFigureNodes->GetNode()); + + m_Controls->m_CircleButton->setChecked(false); + m_Controls->m_RectangleButton->setChecked(false); + m_Controls->m_PolygonButton->setChecked(false); + + //this->GetActiveStdMultiWidget()->SetWidgetPlanesVisibility(true); + this->RequestStatisticsUpdate(); + +} + +void QmitkPartialVolumeAnalysisView::PlanarFigureFocus(mitk::DataNode* node) +{ + mitk::PlanarFigure* _PlanarFigure = 0; + _PlanarFigure = dynamic_cast (node->GetData()); + + if (_PlanarFigure) + { + + FindRenderWindow(node); + + const mitk::PlaneGeometry + * _PlaneGeometry = + dynamic_cast (_PlanarFigure->GetGeometry2D()); + + // make node visible + if (m_SelectedRenderWindow) + { + mitk::Point3D centerP = _PlaneGeometry->GetOrigin(); + m_SelectedRenderWindow->GetSliceNavigationController()->ReorientSlices( + centerP, _PlaneGeometry->GetNormal()); + m_SelectedRenderWindow->GetSliceNavigationController()->SelectSliceByPoint( + centerP); + } + } +} + +void QmitkPartialVolumeAnalysisView::FindRenderWindow(mitk::DataNode* node) +{ + + if(node) + { + mitk::PlanarFigure* _PlanarFigure = 0; + _PlanarFigure = dynamic_cast (node->GetData()); + + if (_PlanarFigure) + { + m_SelectedRenderWindow = 0; + QmitkRenderWindow* RenderWindow1 = + this->GetActiveStdMultiWidget()->GetRenderWindow1(); + QmitkRenderWindow* RenderWindow2 = + this->GetActiveStdMultiWidget()->GetRenderWindow2(); + QmitkRenderWindow* RenderWindow3 = + this->GetActiveStdMultiWidget()->GetRenderWindow3(); + QmitkRenderWindow* RenderWindow4 = + this->GetActiveStdMultiWidget()->GetRenderWindow4(); + + bool PlanarFigureInitializedWindow = false; + + // find initialized renderwindow + if (node->GetBoolProperty("PlanarFigureInitializedWindow", + PlanarFigureInitializedWindow, RenderWindow1->GetRenderer())) + { + m_SelectedRenderWindow = RenderWindow1; + } + + if (!m_SelectedRenderWindow && node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, + RenderWindow2->GetRenderer())) + { + m_SelectedRenderWindow = RenderWindow2; + } + + if (!m_SelectedRenderWindow && node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, + RenderWindow3->GetRenderer())) + { + m_SelectedRenderWindow = RenderWindow3; + } + + if (!m_SelectedRenderWindow && node->GetBoolProperty( + "PlanarFigureInitializedWindow", PlanarFigureInitializedWindow, + RenderWindow4->GetRenderer())) + { + m_SelectedRenderWindow = RenderWindow4; + } + + } + } +} + +void QmitkPartialVolumeAnalysisView::OnSelectionChanged( std::vector nodes ) +{ + + if ( !this->IsVisible() ) + { + return; + } + + if ( nodes.empty() || nodes.size() > 1 ) + { + // Nothing to do: invalidate image, clear statistics, histogram, and GUI + return; + } + + Select(nodes.front()); + +} + +void QmitkPartialVolumeAnalysisView::Select( mitk::DataNode::Pointer node, bool clearMaskOnFirstArgNULL, bool clearImageOnFirstArgNULL ) +{ + // Clear any unreferenced images + this->RemoveOrphanImages(); + + bool somethingChanged = false; + if(node.IsNull()) + { + somethingChanged = true; + + if(clearMaskOnFirstArgNULL) + { + if ( (m_SelectedImageMask.IsNotNull()) && (m_ImageMaskObserverTag >= 0) ) + { + m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); + m_ImageMaskObserverTag = -1; + } + + if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_PlanarFigureObserverTag >= 0) ) + { + m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); + m_PlanarFigureObserverTag = -1; + } + + if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_InitializedObserverTag >= 0) ) + { + m_SelectedPlanarFigure->RemoveObserver( m_InitializedObserverTag ); + m_InitializedObserverTag = -1; + } + + m_SelectedPlanarFigure = NULL; + m_SelectedPlanarFigureNodes->RemoveAllNodes(); + m_CurrentFigureNodeInitialized = false; + m_SelectedRenderWindow = 0; + + m_SelectedMaskNode = NULL; + m_SelectedImageMask = NULL; + + } + + if(clearImageOnFirstArgNULL) + { + if ( (m_SelectedImage.IsNotNull()) && (m_ImageObserverTag >= 0) ) + { + m_SelectedImage->RemoveObserver( m_ImageObserverTag ); + m_ImageObserverTag = -1; + } + + m_SelectedImageNodes->RemoveAllNodes(); + m_SelectedImage = NULL; + + m_IsTensorImage = false; + m_FAImage = NULL; + m_RDImage = NULL; + m_ADImage = NULL; + m_MDImage = NULL; + m_CAImage = NULL; + m_DirectionComp1Image = NULL; + m_DirectionComp2Image = NULL; + m_AngularErrorImage = NULL; + + m_Controls->m_SimilarAnglesFrame->setVisible(false); + m_Controls->m_SimilarAnglesLabel->setVisible(false); + } + } + else + { + typedef itk::SimpleMemberCommand< QmitkPartialVolumeAnalysisView > ITKCommandType; + ITKCommandType::Pointer changeListener; + changeListener = ITKCommandType::New(); + changeListener->SetCallbackFunction( this, &QmitkPartialVolumeAnalysisView::RequestStatisticsUpdate ); + + // Get selected element + mitk::TensorImage *selectedTensorImage = dynamic_cast< mitk::TensorImage * >( node->GetData() ); + mitk::Image *selectedImage = dynamic_cast< mitk::Image * >( node->GetData() ); + mitk::PlanarFigure *selectedPlanar = dynamic_cast< mitk::PlanarFigure * >( node->GetData() ); + + bool isMask = false; + bool isImage = false; + bool isPlanar = false; + bool isTensorImage = false; + + if (selectedTensorImage != NULL) + { + isTensorImage = true; + } + else if(selectedImage != NULL) + { + node->GetPropertyValue("binary", isMask); + isImage = !isMask; + } + else if ( (selectedPlanar != NULL) ) + { + isPlanar = true; + } + + // image + if(isImage && selectedImage->GetDimension()==3) + { + if(selectedImage != m_SelectedImage.GetPointer()) + { + somethingChanged = true; + + if ( (m_SelectedImage.IsNotNull()) && (m_ImageObserverTag >= 0) ) + { + m_SelectedImage->RemoveObserver( m_ImageObserverTag ); + m_ImageObserverTag = -1; + } + + *m_SelectedImageNodes = node; + m_SelectedImage = selectedImage; + m_IsTensorImage = false; + m_FAImage = NULL; + m_RDImage = NULL; + m_ADImage = NULL; + m_MDImage = NULL; + m_CAImage = NULL; + m_DirectionComp1Image = NULL; + m_DirectionComp2Image = NULL; + m_AngularErrorImage = NULL; + + // Add change listeners to selected objects + m_ImageObserverTag = m_SelectedImage->AddObserver( + itk::ModifiedEvent(), changeListener ); + + m_Controls->m_SimilarAnglesFrame->setVisible(false); + m_Controls->m_SimilarAnglesLabel->setVisible(false); + + m_Controls->m_SelectedImageLabel->setText( m_SelectedImageNodes->GetNode()->GetName().c_str() ); + } + } + + //planar + if(isPlanar) + { + if(selectedPlanar != m_SelectedPlanarFigure.GetPointer()) + { + MITK_INFO << "Planar selection changed"; + somethingChanged = true; + + // Possibly previous change listeners + if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_PlanarFigureObserverTag >= 0) ) + { + m_SelectedPlanarFigure->RemoveObserver( m_PlanarFigureObserverTag ); + m_PlanarFigureObserverTag = -1; + } + + if ( (m_SelectedPlanarFigure.IsNotNull()) && (m_InitializedObserverTag >= 0) ) + { + m_SelectedPlanarFigure->RemoveObserver( m_InitializedObserverTag ); + m_InitializedObserverTag = -1; + } + + m_SelectedPlanarFigure = selectedPlanar; + *m_SelectedPlanarFigureNodes = node; + m_CurrentFigureNodeInitialized = selectedPlanar->IsPlaced(); + + m_SelectedMaskNode = NULL; + m_SelectedImageMask = NULL; + + m_PlanarFigureObserverTag = m_SelectedPlanarFigure->AddObserver( + mitk::EndInteractionPlanarFigureEvent(), changeListener ); + + if(!m_CurrentFigureNodeInitialized) + { + typedef itk::SimpleMemberCommand< QmitkPartialVolumeAnalysisView > ITKCommandType; + ITKCommandType::Pointer initializationCommand; + initializationCommand = ITKCommandType::New(); + + // set the callback function of the member command + initializationCommand->SetCallbackFunction( this, &QmitkPartialVolumeAnalysisView::PlanarFigureInitialized ); + + // add an observer + m_InitializedObserverTag = selectedPlanar->AddObserver( mitk::EndPlacementPlanarFigureEvent(), initializationCommand ); + + } + + m_Controls->m_SelectedMaskLabel->setText( m_SelectedPlanarFigureNodes->GetNode()->GetName().c_str() ); + PlanarFigureFocus(node); + } + } + + //mask + if(isMask && selectedImage->GetDimension()==3) + { + if(selectedImage != m_SelectedImage.GetPointer()) + { + somethingChanged = true; + + if ( (m_SelectedImageMask.IsNotNull()) && (m_ImageMaskObserverTag >= 0) ) + { + m_SelectedImageMask->RemoveObserver( m_ImageMaskObserverTag ); + m_ImageMaskObserverTag = -1; + } + + m_SelectedMaskNode = node; + m_SelectedImageMask = selectedImage; + m_SelectedPlanarFigure = NULL; + m_SelectedPlanarFigureNodes->RemoveAllNodes(); + + m_ImageMaskObserverTag = m_SelectedImageMask->AddObserver( + itk::ModifiedEvent(), changeListener ); + + m_Controls->m_SelectedMaskLabel->setText( m_SelectedMaskNode->GetName().c_str() ); + + } + } + + //tensor image + if(isTensorImage && selectedTensorImage->GetDimension()==3) + { + if(selectedImage != m_SelectedImage.GetPointer()) + { + somethingChanged = true; + + if ( (m_SelectedImage.IsNotNull()) && (m_ImageObserverTag >= 0) ) + { + m_SelectedImage->RemoveObserver( m_ImageObserverTag ); + m_ImageObserverTag = -1; + } + + *m_SelectedImageNodes = node; + m_SelectedImage = selectedImage; + + m_IsTensorImage = true; + + ExtractTensorImages(selectedImage); + + // Add change listeners to selected objects + m_ImageObserverTag = m_SelectedImage->AddObserver( + itk::ModifiedEvent(), changeListener ); + + m_Controls->m_SimilarAnglesFrame->setVisible(true); + m_Controls->m_SimilarAnglesLabel->setVisible(true); + + m_Controls->m_SelectedImageLabel->setText( m_SelectedImageNodes->GetNode()->GetName().c_str() ); + } + } + } + + if(somethingChanged) + { + this->SetMeasurementInfoToRenderWindow(""); + + if(m_SelectedPlanarFigure.IsNull() && m_SelectedImageMask.IsNull() ) + { + m_Controls->m_SelectedMaskLabel->setText( "None" ); + m_Controls->m_ResampleOptionsFrame->setEnabled(false); + m_Controls->m_HistogramWidget->setEnabled(false); + m_Controls->m_ClassSelector->setEnabled(false); + m_Controls->m_DisplayHistogramCheckbox->setEnabled(false); + m_Controls->m_AdvancedCheckbox->setEnabled(false); + m_Controls->frame_7->setEnabled(false); + } + else + { + m_Controls->m_ResampleOptionsFrame->setEnabled(true); + m_Controls->m_HistogramWidget->setEnabled(true); + m_Controls->m_ClassSelector->setEnabled(true); + m_Controls->m_DisplayHistogramCheckbox->setEnabled(true); + m_Controls->m_AdvancedCheckbox->setEnabled(true); + m_Controls->frame_7->setEnabled(true); + } + + // Clear statistics / histogram GUI if nothing is selected + if ( m_SelectedImage.IsNull() ) + { + m_Controls->m_PlanarFigureButtonsFrame->setEnabled(false); + m_Controls->m_OpacityFrame->setEnabled(false); + m_Controls->m_SelectedImageLabel->setText( "None" ); + } + else + { + m_Controls->m_PlanarFigureButtonsFrame->setEnabled(true); + m_Controls->m_OpacityFrame->setEnabled(true); + } + + if( m_SelectedImage.IsNull() + || (m_SelectedPlanarFigure.IsNull() && m_SelectedImageMask.IsNull()) ) + { + m_Controls->m_HistogramWidget->ClearItemModel(); + m_CurrentStatisticsValid = false; + m_Controls->m_ErrorMessageLabel->hide(); + } + else + { + this->RequestStatisticsUpdate(); + } + } +} + + +void QmitkPartialVolumeAnalysisView::ShowClusteringResults() +{ + + typedef itk::Image MaskImageType; + mitk::Image::Pointer mask = 0; + MaskImageType::Pointer itkmask = 0; + + if(m_IsTensorImage && m_Controls->m_SimilarAnglesSlider->value() != 0) + { + typedef itk::Image AngularErrorImageType; + typedef mitk::ImageToItk CastType; + CastType::Pointer caster = CastType::New(); + caster->SetInput(m_AngularErrorImage); + caster->Update(); + + typedef itk::BinaryThresholdImageFilter< AngularErrorImageType, MaskImageType > ThreshType; + ThreshType::Pointer thresh = ThreshType::New(); + thresh->SetUpperThreshold((90-m_Controls->m_SimilarAnglesSlider->value())*(PVA_PI/180.0)); + thresh->SetInsideValue(1.0); + thresh->SetInput(caster->GetOutput()); + thresh->Update(); + itkmask = thresh->GetOutput(); + + mask = mitk::Image::New(); + mask->InitializeByItk(itkmask.GetPointer()); + mask->SetVolume(itkmask->GetBufferPointer()); + + // GetDefaultDataStorage()->Remove(m_newnode); + // m_newnode = mitk::DataNode::New(); + // m_newnode->SetData(mask); + // m_newnode->SetName("masking node"); + // m_newnode->SetIntProperty( "layer", 1002 ); + // GetDefaultDataStorage()->Add(m_newnode, m_SelectedImageNodes->GetNode()); + + } + + mitk::Image::Pointer clusteredImage; + ClusteringType::Pointer clusterer = ClusteringType::New(); + + if(m_QuantifyClass==3) + { + + if(m_IsTensorImage) + { + double *green_fa, *green_rd, *green_ad, *green_md; + //double *greengray_fa, *greengray_rd, *greengray_ad, *greengray_md; + double *gray_fa, *gray_rd, *gray_ad, *gray_md; + //double *redgray_fa, *redgray_rd, *redgray_ad, *redgray_md; + double *red_fa, *red_rd, *red_ad, *red_md; + + mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(0); + mitk::Image::ConstPointer imgToCluster = tmpImg; + + red_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->r, mask); + green_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->g, mask); + gray_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->b, mask); + + tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(3); + mitk::Image::ConstPointer imgToCluster3 = tmpImg; + + red_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentRGBClusteringResults->rgbChannels->r, mask); + green_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentRGBClusteringResults->rgbChannels->g, mask); + gray_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentRGBClusteringResults->rgbChannels->b, mask); + + tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(4); + mitk::Image::ConstPointer imgToCluster4 = tmpImg; + + red_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentRGBClusteringResults->rgbChannels->r, mask); + green_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentRGBClusteringResults->rgbChannels->g, mask); + gray_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentRGBClusteringResults->rgbChannels->b, mask); + + tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(5); + mitk::Image::ConstPointer imgToCluster5 = tmpImg; + + red_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentRGBClusteringResults->rgbChannels->r, mask); + green_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentRGBClusteringResults->rgbChannels->g, mask); + gray_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentRGBClusteringResults->rgbChannels->b, mask); + + // clipboard + QString clipboardText("FA\t%1\t%2\t\t%3\t%4\t\t%5\t%6\t"); + clipboardText = clipboardText + .arg(red_fa[0]).arg(red_fa[1]) + .arg(gray_fa[0]).arg(gray_fa[1]) + .arg(green_fa[0]).arg(green_fa[1]); + QString clipboardText3("RD\t%1\t%2\t\t%3\t%4\t\t%5\t%6\t"); + clipboardText3 = clipboardText3 + .arg(red_rd[0]).arg(red_rd[1]) + .arg(gray_rd[0]).arg(gray_rd[1]) + .arg(green_rd[0]).arg(green_rd[1]); + QString clipboardText4("AD\t%1\t%2\t\t%3\t%4\t\t%5\t%6\t"); + clipboardText4 = clipboardText4 + .arg(red_ad[0]).arg(red_ad[1]) + .arg(gray_ad[0]).arg(gray_ad[1]) + .arg(green_ad[0]).arg(green_ad[1]); + QString clipboardText5("MD\t%1\t%2\t\t%3\t%4\t\t%5\t%6"); + clipboardText5 = clipboardText5 + .arg(red_md[0]).arg(red_md[1]) + .arg(gray_md[0]).arg(gray_md[1]) + .arg(green_md[0]).arg(green_md[1]); + + QApplication::clipboard()->setText(clipboardText+clipboardText3+clipboardText4+clipboardText5, QClipboard::Clipboard); + + // now paint infos also on renderwindow + QString plainInfoText("%1 %2 %3 \n"); + plainInfoText = plainInfoText + .arg("Red ", 20) + .arg("Gray ", 20) + .arg("Green", 20); + + QString plainInfoText0("FA:%1 ± %2%3 ± %4%5 ± %6\n"); + plainInfoText0 = plainInfoText0 + .arg(red_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(red_fa[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(gray_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(gray_fa[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(green_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(green_fa[1], -10, 'g', 2, QLatin1Char( ' ' )); + + QString plainInfoText3("RDx10³:%1 ± %2%3 ± %4%5 ± %6\n"); + plainInfoText3 = plainInfoText3 + .arg(1000.0 * red_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_rd[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(1000.0 * gray_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * gray_rd[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(1000.0 * green_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * green_rd[1], -10, 'g', 2, QLatin1Char( ' ' )); + + QString plainInfoText4("ADx10³:%1 ± %2%3 ± %4%5 ± %6\n"); + plainInfoText4 = plainInfoText4 + .arg(1000.0 * red_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_ad[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(1000.0 * gray_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * gray_ad[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(1000.0 * green_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * green_ad[1], -10, 'g', 2, QLatin1Char( ' ' )); + + QString plainInfoText5("MDx10³:%1 ± %2%3 ± %4%5 ± %6"); + plainInfoText5 = plainInfoText5 + .arg(1000.0 * red_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_md[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(1000.0 * gray_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * gray_md[1], -10, 'g', 2, QLatin1Char( ' ' )) + .arg(1000.0 * green_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * green_md[1], -10, 'g', 2, QLatin1Char( ' ' )); + + this->SetMeasurementInfoToRenderWindow(plainInfoText+plainInfoText0+plainInfoText3+plainInfoText4+plainInfoText5); + + } + else + { + double* green; + double* gray; + double* red; + + mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalImage(); + mitk::Image::ConstPointer imgToCluster = tmpImg; + + red = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->r); + green = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->g); + gray = clusterer->PerformQuantification(imgToCluster, m_CurrentRGBClusteringResults->rgbChannels->b); + + // clipboard + QString clipboardText("%1\t%2\t\t%3\t%4\t\t%5\t%6"); + clipboardText = clipboardText.arg(red[0]).arg(red[1]) + .arg(gray[0]).arg(gray[1]) + .arg(green[0]).arg(green[1]); + QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); + + // now paint infos also on renderwindow + QString plainInfoText("Red: %1 ± %2\nGray: %3 ± %4\nGreen: %5 ± %6"); + plainInfoText = plainInfoText.arg(red[0]).arg(red[1]) + .arg(gray[0]).arg(gray[1]) + .arg(green[0]).arg(green[1]); + + this->SetMeasurementInfoToRenderWindow(plainInfoText); + + } + + + clusteredImage = m_CurrentRGBClusteringResults->rgb; + + } + else + { + if(m_IsTensorImage) + { + double *red_fa, *red_rd, *red_ad, *red_md; + + mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(0); + mitk::Image::ConstPointer imgToCluster = tmpImg; + red_fa = clusterer->PerformQuantification(imgToCluster, m_CurrentPerformClusteringResults->clusteredImage, mask); + + tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(3); + mitk::Image::ConstPointer imgToCluster3 = tmpImg; + red_rd = clusterer->PerformQuantification(imgToCluster3, m_CurrentPerformClusteringResults->clusteredImage, mask); + + tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(4); + mitk::Image::ConstPointer imgToCluster4 = tmpImg; + red_ad = clusterer->PerformQuantification(imgToCluster4, m_CurrentPerformClusteringResults->clusteredImage, mask); + + tmpImg = m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(5); + mitk::Image::ConstPointer imgToCluster5 = tmpImg; + red_md = clusterer->PerformQuantification(imgToCluster5, m_CurrentPerformClusteringResults->clusteredImage, mask); + + // clipboard + QString clipboardText("FA\t%1\t%2\t"); + clipboardText = clipboardText + .arg(red_fa[0]).arg(red_fa[1]); + QString clipboardText3("RD\t%1\t%2\t"); + clipboardText3 = clipboardText3 + .arg(red_rd[0]).arg(red_rd[1]); + QString clipboardText4("AD\t%1\t%2\t"); + clipboardText4 = clipboardText4 + .arg(red_ad[0]).arg(red_ad[1]); + QString clipboardText5("MD\t%1\t%2\t"); + clipboardText5 = clipboardText5 + .arg(red_md[0]).arg(red_md[1]); + + QApplication::clipboard()->setText(clipboardText+clipboardText3+clipboardText4+clipboardText5, QClipboard::Clipboard); + + // now paint infos also on renderwindow + QString plainInfoText("%1 \n"); + plainInfoText = plainInfoText + .arg("Red ", 20); + + QString plainInfoText0("FA:%1 ± %2\n"); + plainInfoText0 = plainInfoText0 + .arg(red_fa[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(red_fa[1], -10, 'g', 2, QLatin1Char( ' ' )); + + QString plainInfoText3("RDx10³:%1 ± %2\n"); + plainInfoText3 = plainInfoText3 + .arg(1000.0 * red_rd[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_rd[1], -10, 'g', 2, QLatin1Char( ' ' )); + + QString plainInfoText4("ADx10³:%1 ± %2\n"); + plainInfoText4 = plainInfoText4 + .arg(1000.0 * red_ad[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_ad[1], -10, 'g', 2, QLatin1Char( ' ' )); + + QString plainInfoText5("MDx10³:%1 ± %2"); + plainInfoText5 = plainInfoText5 + .arg(1000.0 * red_md[0], 10, 'g', 2, QLatin1Char( ' ' )).arg(1000.0 * red_md[1], -10, 'g', 2, QLatin1Char( ' ' )); + + this->SetMeasurementInfoToRenderWindow(plainInfoText+plainInfoText0+plainInfoText3+plainInfoText4+plainInfoText5); + + } + else + { + double* quant; + + mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalImage(); + mitk::Image::ConstPointer imgToCluster = tmpImg; + quant = clusterer->PerformQuantification(imgToCluster, m_CurrentPerformClusteringResults->clusteredImage); + + // clipboard + QString clipboardText("%1\t%2"); + clipboardText = clipboardText.arg(quant[0]).arg(quant[1]); + QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); + + // now paint infos also on renderwindow + QString plainInfoText("Measurement: %1 ± %2"); + plainInfoText = plainInfoText.arg(quant[0]).arg(quant[1]); + + this->SetMeasurementInfoToRenderWindow(plainInfoText); + } + + clusteredImage = m_CurrentPerformClusteringResults->displayImage; + + } + + if(mask.IsNotNull()) + { + typedef itk::Image,3> RGBImageType; + typedef mitk::ImageToItk ClusterCasterType; + ClusterCasterType::Pointer clCaster = ClusterCasterType::New(); + clCaster->SetInput(clusteredImage); + clCaster->Update(); + clCaster->GetOutput(); + + typedef itk::MaskImageFilter< RGBImageType, MaskImageType, RGBImageType > MaskType; + MaskType::Pointer masker = MaskType::New(); + masker->SetInput1(clCaster->GetOutput()); + masker->SetInput2(itkmask); + masker->Update(); + + clusteredImage = mitk::Image::New(); + clusteredImage->InitializeByItk(masker->GetOutput()); + clusteredImage->SetVolume(masker->GetOutput()->GetBufferPointer()); + } + + if(m_ClusteringResult.IsNotNull()) + { + GetDefaultDataStorage()->Remove(m_ClusteringResult); + } + + m_ClusteringResult = mitk::DataNode::New(); + m_ClusteringResult->SetBoolProperty("helper object", true); + m_ClusteringResult->SetIntProperty( "layer", 1000 ); + m_ClusteringResult->SetBoolProperty("texture interpolation", m_TexIsOn); + m_ClusteringResult->SetData(clusteredImage); + m_ClusteringResult->SetName("Clusterprobs"); + GetDefaultDataStorage()->Add(m_ClusteringResult, m_SelectedImageNodes->GetNode()); + + if(m_SelectedPlanarFigure.IsNotNull() && m_SelectedPlanarFigureNodes->GetNode().IsNotNull()) + { + m_SelectedPlanarFigureNodes->GetNode()->SetIntProperty( "layer", 1001 ); + } + + GetActiveStdMultiWidget()->RequestUpdate(); + +} + + +void QmitkPartialVolumeAnalysisView::UpdateStatistics() +{ + + MITK_INFO << "UpdateStatistics()"; + + if(!m_CurrentFigureNodeInitialized && m_SelectedPlanarFigure.IsNotNull()) + { + MITK_INFO << "Selected planar figure not initialized. No stats calculation performed."; + return; + } + + // Remove any cached images that are no longer referenced elsewhere + this->RemoveOrphanImages(); + + QmitkStdMultiWidget *multiWidget = this->GetActiveStdMultiWidget(); + if ( multiWidget == NULL ) + { + return; + } + + if ( m_SelectedImage.IsNotNull() ) + { + // Check if a the selected image is a multi-channel image. If yes, statistics + // cannot be calculated currently. + if ( !m_IsTensorImage && m_SelectedImage->GetPixelType().GetNumberOfComponents() > 1 ) + { + std::stringstream message; + message << "Multi-component images not supported."; + m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); + m_Controls->m_ErrorMessageLabel->show(); + + m_Controls->m_HistogramWidget->ClearItemModel(); + m_CurrentStatisticsValid = false; + return; + } + + // Retrieve HistogramStatisticsCalculator from has map (or create a new one + // for this image if non-existant) + PartialVolumeAnalysisMapType::iterator it = + m_PartialVolumeAnalysisMap.find( m_SelectedImage ); + + if ( it != m_PartialVolumeAnalysisMap.end() ) + { + m_CurrentStatisticsCalculator = it->second; + MITK_INFO << "Retrieving StatisticsCalculator"; + } + else + { + m_CurrentStatisticsCalculator = mitk::PartialVolumeAnalysisHistogramCalculator::New(); + if(m_IsTensorImage) + { + m_CurrentStatisticsCalculator->SetImage( m_CAImage ); + m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_FAImage ); + m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_DirectionComp1Image ); + m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_DirectionComp2Image ); + m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_RDImage ); + m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_ADImage ); + m_CurrentStatisticsCalculator->AddAdditionalResamplingImage( m_MDImage ); + } + else + { + m_CurrentStatisticsCalculator->SetImage( m_SelectedImage ); + } + m_PartialVolumeAnalysisMap[m_SelectedImage] = m_CurrentStatisticsCalculator; + MITK_INFO << "Creating StatisticsCalculator"; + } + + std::string maskName; + std::string maskType; + unsigned int maskDimension; + + if ( m_SelectedImageMask.IsNotNull() ) + { + m_CurrentStatisticsCalculator->SetImageMask( m_SelectedImageMask ); + m_CurrentStatisticsCalculator->SetMaskingModeToImage(); + + maskName = m_SelectedMaskNode->GetName(); + maskType = m_SelectedImageMask->GetNameOfClass(); + maskDimension = 3; + + std::stringstream maskLabel; + maskLabel << maskName; + if ( maskDimension > 0 ) + { + maskLabel << " [" << maskDimension << "D " << maskType << "]"; + } + m_Controls->m_SelectedMaskLabel->setText( maskLabel.str().c_str() ); + } + else if ( m_SelectedPlanarFigure.IsNotNull() && m_SelectedPlanarFigureNodes->GetNode().IsNotNull()) + { + m_CurrentStatisticsCalculator->SetPlanarFigure( m_SelectedPlanarFigure ); + m_CurrentStatisticsCalculator->SetMaskingModeToPlanarFigure(); + + maskName = m_SelectedPlanarFigureNodes->GetNode()->GetName(); + maskType = m_SelectedPlanarFigure->GetNameOfClass(); + maskDimension = 2; + } + else + { + m_CurrentStatisticsCalculator->SetMaskingModeToNone(); + + maskName = "None"; + maskType = ""; + maskDimension = 0; + } + + bool statisticsChanged = false; + bool statisticsCalculationSuccessful = false; + + // Initialize progress bar + mitk::ProgressBar::GetInstance()->AddStepsToDo( 100 ); + + // Install listener for progress events and initialize progress bar + typedef itk::SimpleMemberCommand< QmitkPartialVolumeAnalysisView > ITKCommandType; + ITKCommandType::Pointer progressListener; + progressListener = ITKCommandType::New(); + progressListener->SetCallbackFunction( this, &QmitkPartialVolumeAnalysisView::UpdateProgressBar ); + unsigned long progressObserverTag = m_CurrentStatisticsCalculator + ->AddObserver( itk::ProgressEvent(), progressListener ); + + ClusteringType::ParamsType *cparams = 0; + ClusteringType::ClusterResultType *cresult = 0; + ClusteringType::HistType *chist = 0; + + try + { + m_CurrentStatisticsCalculator->SetNumberOfBins(m_Controls->m_NumberBins->text().toInt()); + m_CurrentStatisticsCalculator->SetUpsamplingFactor(m_Controls->m_Upsampling->text().toDouble()); + m_CurrentStatisticsCalculator->SetGaussianSigma(m_Controls->m_GaussianSigma->text().toDouble()); + + // Compute statistics + statisticsChanged = + m_CurrentStatisticsCalculator->ComputeStatistics( ); + + mitk::Image* tmpImg = m_CurrentStatisticsCalculator->GetInternalImage(); + mitk::Image::ConstPointer imgToCluster = tmpImg; + if(imgToCluster.IsNotNull()) + { + + // perform clustering + const HistogramType *histogram = m_CurrentStatisticsCalculator->GetHistogram( ); + ClusteringType::Pointer clusterer = ClusteringType::New(); + clusterer->SetStepsNumIntegration(200); + clusterer->SetMaxIt(1000); + + mitk::Image::Pointer pFiberImg; + if(m_QuantifyClass==3) + { + if(m_Controls->m_Quantiles->isChecked()) + { + m_CurrentRGBClusteringResults = clusterer->PerformRGBQuantiles(imgToCluster, histogram, m_Controls->m_q1->value(),m_Controls->m_q2->value()); + } + else + { + m_CurrentRGBClusteringResults = clusterer->PerformRGBClustering(imgToCluster, histogram); + } + + pFiberImg = m_CurrentRGBClusteringResults->rgbChannels->r; + cparams = m_CurrentRGBClusteringResults->params; + cresult = m_CurrentRGBClusteringResults->result; + chist = m_CurrentRGBClusteringResults->hist; + } + else + { + if(m_Controls->m_Quantiles->isChecked()) + { + m_CurrentPerformClusteringResults = + clusterer->PerformQuantiles(imgToCluster, histogram, m_Controls->m_q1->value(),m_Controls->m_q2->value()); + } + else + { + m_CurrentPerformClusteringResults = + clusterer->PerformClustering(imgToCluster, histogram, m_QuantifyClass); + } + + pFiberImg = m_CurrentPerformClusteringResults->clusteredImage; + cparams = m_CurrentPerformClusteringResults->params; + cresult = m_CurrentPerformClusteringResults->result; + chist = m_CurrentPerformClusteringResults->hist; + } + + if(m_IsTensorImage) + { + m_AngularErrorImage = clusterer->CaculateAngularErrorImage( + m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(1), + m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(2), + pFiberImg); + + // GetDefaultDataStorage()->Remove(m_newnode2); + // m_newnode2 = mitk::DataNode::New(); + // m_newnode2->SetData(m_AngularErrorImage); + // m_newnode2->SetName(("AngularError")); + // m_newnode2->SetIntProperty( "layer", 1003 ); + // GetDefaultDataStorage()->Add(m_newnode2, m_SelectedImageNodes->GetNode()); + + // newnode = mitk::DataNode::New(); + // newnode->SetData(m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(1)); + // newnode->SetName(("Comp1")); + // GetDefaultDataStorage()->Add(newnode, m_SelectedImageNodes->GetNode()); + + // newnode = mitk::DataNode::New(); + // newnode->SetData(m_CurrentStatisticsCalculator->GetInternalAdditionalResampledImage(2)); + // newnode->SetName(("Comp2")); + // GetDefaultDataStorage()->Add(newnode, m_SelectedImageNodes->GetNode()); + } + ShowClusteringResults(); + } + + statisticsCalculationSuccessful = true; + } + catch ( const std::runtime_error &e ) + { + // In case of exception, print error message on GUI + std::stringstream message; + message << "" << e.what() << ""; + m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); + m_Controls->m_ErrorMessageLabel->show(); + } + catch ( const std::exception &e ) + { + MITK_ERROR << "Caught exception: " << e.what(); + + // In case of exception, print error message on GUI + std::stringstream message; + message << "Error in calculating histogram: " << e.what() << ""; + m_Controls->m_ErrorMessageLabel->setText( message.str().c_str() ); + m_Controls->m_ErrorMessageLabel->show(); + } + + m_CurrentStatisticsCalculator->RemoveObserver( progressObserverTag ); + + // Make sure that progress bar closes + mitk::ProgressBar::GetInstance()->Progress( 100 ); + + if ( statisticsCalculationSuccessful ) + { + if ( statisticsChanged ) + { + // Do not show any error messages + m_Controls->m_ErrorMessageLabel->hide(); + + m_CurrentStatisticsValid = true; + } + + // m_Controls->m_HistogramWidget->SetHistogramModeToDirectHistogram(); + m_Controls->m_HistogramWidget->SetParameters( + cparams, cresult, chist ); + // m_Controls->m_HistogramWidget->UpdateItemModelFromHistogram(); + + } + else + { + m_Controls->m_SelectedMaskLabel->setText( "None" ); + + // Clear statistics and histogram + m_Controls->m_HistogramWidget->ClearItemModel(); + m_CurrentStatisticsValid = false; + + + // If a (non-closed) PlanarFigure is selected, display a line profile widget + if ( m_SelectedPlanarFigure.IsNotNull() ) + { + // TODO: enable line profile widget + //m_Controls->m_StatisticsWidgetStack->setCurrentIndex( 1 ); + //m_Controls->m_LineProfileWidget->SetImage( m_SelectedImage ); + //m_Controls->m_LineProfileWidget->SetPlanarFigure( m_SelectedPlanarFigure ); + //m_Controls->m_LineProfileWidget->UpdateItemModelFromPath(); + } + } + } + +} + +void QmitkPartialVolumeAnalysisView::SetMeasurementInfoToRenderWindow(const QString& text) +{ + + FindRenderWindow(m_SelectedPlanarFigureNodes->GetNode()); + + if(m_LastRenderWindow != m_SelectedRenderWindow) + { + + if(m_LastRenderWindow) + { + QObject::disconnect( m_LastRenderWindow, SIGNAL( destroyed(QObject*) ) + , this, SLOT( OnRenderWindowDelete(QObject*) ) ); + } + m_LastRenderWindow = m_SelectedRenderWindow; + if(m_LastRenderWindow) + { + QObject::connect( m_LastRenderWindow, SIGNAL( destroyed(QObject*) ) + , this, SLOT( OnRenderWindowDelete(QObject*) ) ); + } + } + + if(m_LastRenderWindow && m_SelectedPlanarFigureNodes->GetNode().IsNotNull()) + { + if (!text.isEmpty()) + { + m_MeasurementInfoAnnotation->SetText(1, text.toLatin1().data()); + mitk::VtkLayerController::GetInstance(m_LastRenderWindow->GetRenderWindow())->InsertForegroundRenderer( + m_MeasurementInfoRenderer, true); + } + else + { + if (mitk::VtkLayerController::GetInstance( + m_LastRenderWindow->GetRenderWindow()) ->IsRendererInserted( + m_MeasurementInfoRenderer)) + mitk::VtkLayerController::GetInstance(m_LastRenderWindow->GetRenderWindow())->RemoveRenderer( + m_MeasurementInfoRenderer); + } + } + else + { + if (!text.isEmpty()) + { + m_MeasurementInfoAnnotation->SetText(1, text.toLatin1().data()); + mitk::VtkLayerController::GetInstance(this->GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderWindow())->InsertForegroundRenderer( + m_MeasurementInfoRenderer, true); + } + else + { + if (mitk::VtkLayerController::GetInstance( + this->GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderWindow()) ->IsRendererInserted( + m_MeasurementInfoRenderer)) + mitk::VtkLayerController::GetInstance(this->GetActiveStdMultiWidget()->GetRenderWindow1()->GetRenderWindow())->RemoveRenderer( + m_MeasurementInfoRenderer); + } + } +} + +void QmitkPartialVolumeAnalysisView::UpdateProgressBar() +{ + mitk::ProgressBar::GetInstance()->Progress(); +} + +void QmitkPartialVolumeAnalysisView::RequestStatisticsUpdate() +{ + if ( !m_StatisticsUpdatePending ) + { + QApplication::postEvent( this, new QmitkRequestStatisticsUpdateEvent ); + m_StatisticsUpdatePending = true; + } +} + + +void QmitkPartialVolumeAnalysisView::RemoveOrphanImages() +{ + PartialVolumeAnalysisMapType::iterator it = m_PartialVolumeAnalysisMap.begin(); + + while ( it != m_PartialVolumeAnalysisMap.end() ) + { + mitk::Image *image = it->first; + mitk::PartialVolumeAnalysisHistogramCalculator *calculator = it->second; + ++it; + + mitk::NodePredicateData::Pointer hasImage = mitk::NodePredicateData::New( image ); + if ( this->GetDefaultDataStorage()->GetNode( hasImage ) == NULL ) + { + if ( m_SelectedImage == image ) + { + m_SelectedImage = NULL; + m_SelectedImageNodes->RemoveAllNodes(); + } + if ( m_CurrentStatisticsCalculator == calculator ) + { + m_CurrentStatisticsCalculator = NULL; + } + m_PartialVolumeAnalysisMap.erase( image ); + it = m_PartialVolumeAnalysisMap.begin(); + } + } +} + +void QmitkPartialVolumeAnalysisView::ExtractTensorImages( + mitk::Image::ConstPointer tensorimage) +{ + typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; + + typedef mitk::ImageToItk CastType; + CastType::Pointer caster = CastType::New(); + caster->SetInput(tensorimage); + caster->Update(); + TensorImageType::Pointer image = caster->GetOutput(); + + typedef itk::TensorDerivedMeasurementsFilter MeasurementsType; + MeasurementsType::Pointer measurementsCalculator = MeasurementsType::New(); + measurementsCalculator->SetInput(image ); + measurementsCalculator->SetMeasure(MeasurementsType::FA); + measurementsCalculator->Update(); + MeasurementsType::OutputImageType::Pointer fa = measurementsCalculator->GetOutput(); + + m_FAImage = mitk::Image::New(); + m_FAImage->InitializeByItk(fa.GetPointer()); + m_FAImage->SetVolume(fa->GetBufferPointer()); + + // mitk::DataNode::Pointer node = mitk::DataNode::New(); + // node->SetData(m_FAImage); + // GetDefaultDataStorage()->Add(node); + + measurementsCalculator = MeasurementsType::New(); + measurementsCalculator->SetInput(image ); + measurementsCalculator->SetMeasure(MeasurementsType::CA); + measurementsCalculator->Update(); + MeasurementsType::OutputImageType::Pointer ca = measurementsCalculator->GetOutput(); + + m_CAImage = mitk::Image::New(); + m_CAImage->InitializeByItk(ca.GetPointer()); + m_CAImage->SetVolume(ca->GetBufferPointer()); + + // node = mitk::DataNode::New(); + // node->SetData(m_CAImage); + // GetDefaultDataStorage()->Add(node); + + measurementsCalculator = MeasurementsType::New(); + measurementsCalculator->SetInput(image ); + measurementsCalculator->SetMeasure(MeasurementsType::RD); + measurementsCalculator->Update(); + MeasurementsType::OutputImageType::Pointer rd = measurementsCalculator->GetOutput(); + + m_RDImage = mitk::Image::New(); + m_RDImage->InitializeByItk(rd.GetPointer()); + m_RDImage->SetVolume(rd->GetBufferPointer()); + + // node = mitk::DataNode::New(); + // node->SetData(m_CAImage); + // GetDefaultDataStorage()->Add(node); + + measurementsCalculator = MeasurementsType::New(); + measurementsCalculator->SetInput(image ); + measurementsCalculator->SetMeasure(MeasurementsType::AD); + measurementsCalculator->Update(); + MeasurementsType::OutputImageType::Pointer ad = measurementsCalculator->GetOutput(); + + m_ADImage = mitk::Image::New(); + m_ADImage->InitializeByItk(ad.GetPointer()); + m_ADImage->SetVolume(ad->GetBufferPointer()); + + // node = mitk::DataNode::New(); + // node->SetData(m_CAImage); + // GetDefaultDataStorage()->Add(node); + + measurementsCalculator = MeasurementsType::New(); + measurementsCalculator->SetInput(image ); + measurementsCalculator->SetMeasure(MeasurementsType::RA); + measurementsCalculator->Update(); + MeasurementsType::OutputImageType::Pointer md = measurementsCalculator->GetOutput(); + + m_MDImage = mitk::Image::New(); + m_MDImage->InitializeByItk(md.GetPointer()); + m_MDImage->SetVolume(md->GetBufferPointer()); + + // node = mitk::DataNode::New(); + // node->SetData(m_CAImage); + // GetDefaultDataStorage()->Add(node); + + typedef DirectionsFilterType::OutputImageType DirImageType; + DirectionsFilterType::Pointer dirFilter = DirectionsFilterType::New(); + dirFilter->SetInput(image ); + dirFilter->Update(); + + itk::ImageRegionIterator + itd(dirFilter->GetOutput(), dirFilter->GetOutput()->GetLargestPossibleRegion()); + itd = itd.Begin(); + while( !itd.IsAtEnd() ) + { + DirImageType::PixelType direction = itd.Get(); + direction[0] = fabs(direction[0]); + direction[1] = fabs(direction[1]); + direction[2] = fabs(direction[2]); + itd.Set(direction); + ++itd; + } + + typedef itk::CartesianToPolarVectorImageFilter< + DirImageType, DirImageType, true> C2PFilterType; + C2PFilterType::Pointer cpFilter = C2PFilterType::New(); + cpFilter->SetInput(dirFilter->GetOutput()); + cpFilter->Update(); + DirImageType::Pointer dir = cpFilter->GetOutput(); + + typedef itk::Image CompImageType; + CompImageType::Pointer comp1 = CompImageType::New(); + comp1->SetSpacing( dir->GetSpacing() ); // Set the image spacing + comp1->SetOrigin( dir->GetOrigin() ); // Set the image origin + comp1->SetDirection( dir->GetDirection() ); // Set the image direction + comp1->SetRegions( dir->GetLargestPossibleRegion() ); + comp1->Allocate(); + + CompImageType::Pointer comp2 = CompImageType::New(); + comp2->SetSpacing( dir->GetSpacing() ); // Set the image spacing + comp2->SetOrigin( dir->GetOrigin() ); // Set the image origin + comp2->SetDirection( dir->GetDirection() ); // Set the image direction + comp2->SetRegions( dir->GetLargestPossibleRegion() ); + comp2->Allocate(); + + itk::ImageRegionConstIterator + it(dir, dir->GetLargestPossibleRegion()); + + itk::ImageRegionIterator + it1(comp1, comp1->GetLargestPossibleRegion()); + + itk::ImageRegionIterator + it2(comp2, comp2->GetLargestPossibleRegion()); + + it = it.Begin(); + it1 = it1.Begin(); + it2 = it2.Begin(); + + while( !it.IsAtEnd() ) + { + it1.Set(it.Get()[1]); + it2.Set(it.Get()[2]); + ++it; + ++it1; + ++it2; + } + + m_DirectionComp1Image = mitk::Image::New(); + m_DirectionComp1Image->InitializeByItk(comp1.GetPointer()); + m_DirectionComp1Image->SetVolume(comp1->GetBufferPointer()); + + m_DirectionComp2Image = mitk::Image::New(); + m_DirectionComp2Image->InitializeByItk(comp2.GetPointer()); + m_DirectionComp2Image->SetVolume(comp2->GetBufferPointer()); + +} + +void QmitkPartialVolumeAnalysisView::OnRenderWindowDelete(QObject * obj) +{ + if(obj == m_LastRenderWindow) + m_LastRenderWindow = 0; +} + +bool QmitkPartialVolumeAnalysisView::event( QEvent *event ) +{ + if ( event->type() == (QEvent::Type) QmitkRequestStatisticsUpdateEvent::StatisticsUpdateRequest ) + { + // Update statistics + + m_StatisticsUpdatePending = false; + + this->UpdateStatistics(); + return true; + } + + return false; +} + +void QmitkPartialVolumeAnalysisView::Visible() +{ + this->OnSelectionChanged( this->GetDataManagerSelection() ); +} + +bool QmitkPartialVolumeAnalysisView::IsExclusiveFunctionality() const +{ + return true; +} + +void QmitkPartialVolumeAnalysisView::Activated() +{ + this->GetActiveStdMultiWidget()->SetWidgetPlanesVisibility(false); + //this->GetActiveStdMultiWidget()->GetRenderWindow1()->FullScreenMode(true); + + mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); + mitk::DataNode* node = 0; + mitk::PlanarFigure* figure = 0; + mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; + + // finally add all nodes to the model + for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() + ; it++) + { + node = const_cast(it->Value().GetPointer()); + figure = dynamic_cast(node->GetData()); + + if(figure) + { + figureInteractor = dynamic_cast(node->GetInteractor()); + + if(figureInteractor.IsNull()) + figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", node); + + mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); + } + } + + m_Visible = true; + +} + +void QmitkPartialVolumeAnalysisView::Deactivated() +{ + this->GetActiveStdMultiWidget()->SetWidgetPlanesVisibility(true); + //this->GetActiveStdMultiWidget()->GetRenderWindow1()->FullScreenMode(false); + + this->SetMeasurementInfoToRenderWindow(""); + + mitk::DataStorage::SetOfObjects::ConstPointer _NodeSet = this->GetDefaultDataStorage()->GetAll(); + mitk::DataNode* node = 0; + mitk::PlanarFigure* figure = 0; + mitk::PlanarFigureInteractor::Pointer figureInteractor = 0; + + // finally add all nodes to the model + for(mitk::DataStorage::SetOfObjects::ConstIterator it=_NodeSet->Begin(); it!=_NodeSet->End() + ; it++) + { + node = const_cast(it->Value().GetPointer()); + figure = dynamic_cast(node->GetData()); + + if(figure) + { + figureInteractor = dynamic_cast(node->GetInteractor()); + + if(figureInteractor) + mitk::GlobalInteraction::GetInstance()->RemoveInteractor(figureInteractor); + } + } + + m_Visible = false; +} + +void QmitkPartialVolumeAnalysisView::GreenRadio(bool checked) +{ + if(checked) + { + m_Controls->m_PartialVolumeRadio->setChecked(false); + m_Controls->m_BlueRadio->setChecked(false); + m_Controls->m_AllRadio->setChecked(false); + } + + m_QuantifyClass = 0; + + RequestStatisticsUpdate(); +} + + +void QmitkPartialVolumeAnalysisView::PartialVolumeRadio(bool checked) +{ + if(checked) + { + m_Controls->m_GreenRadio->setChecked(false); + m_Controls->m_BlueRadio->setChecked(false); + m_Controls->m_AllRadio->setChecked(false); + } + + m_QuantifyClass = 1; + + RequestStatisticsUpdate(); +} + +void QmitkPartialVolumeAnalysisView::BlueRadio(bool checked) +{ + if(checked) + { + m_Controls->m_PartialVolumeRadio->setChecked(false); + m_Controls->m_GreenRadio->setChecked(false); + m_Controls->m_AllRadio->setChecked(false); + } + + m_QuantifyClass = 2; + + RequestStatisticsUpdate(); +} + +void QmitkPartialVolumeAnalysisView::AllRadio(bool checked) +{ + if(checked) + { + m_Controls->m_BlueRadio->setChecked(false); + m_Controls->m_PartialVolumeRadio->setChecked(false); + m_Controls->m_GreenRadio->setChecked(false); + } + + m_QuantifyClass = 3; + + RequestStatisticsUpdate(); +} + +void QmitkPartialVolumeAnalysisView::NumberBinsChangedSlider(int v ) +{ + m_Controls->m_NumberBins->setText(QString("%1").arg(m_Controls->m_NumberBinsSlider->value()*5.0)); +} + +void QmitkPartialVolumeAnalysisView::UpsamplingChangedSlider( int v) +{ + m_Controls->m_Upsampling->setText(QString("%1").arg(m_Controls->m_UpsamplingSlider->value()/10.0)); +} + +void QmitkPartialVolumeAnalysisView::GaussianSigmaChangedSlider(int v ) +{ + m_Controls->m_GaussianSigma->setText(QString("%1").arg(m_Controls->m_GaussianSigmaSlider->value()/100.0)); +} + +void QmitkPartialVolumeAnalysisView::SimilarAnglesChangedSlider(int v ) +{ + m_Controls->m_SimilarAngles->setText(QString("%1°").arg(90-m_Controls->m_SimilarAnglesSlider->value())); + ShowClusteringResults(); +} + +void QmitkPartialVolumeAnalysisView::OpacityChangedSlider(int v ) +{ + + if(m_SelectedImageNodes->GetNode().IsNotNull()) + { + float opacImag = 1.0f-(v-5)/5.0f; + opacImag = opacImag < 0 ? 0 : opacImag; + m_SelectedImageNodes->GetNode()->SetFloatProperty("opacity", opacImag); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } + + if(m_ClusteringResult.IsNotNull()) + { + float opacClust = v/5.0f; + opacClust = opacClust > 1 ? 1 : opacClust; + m_ClusteringResult->SetFloatProperty("opacity", opacClust); + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); + } +} + +void QmitkPartialVolumeAnalysisView::NumberBinsReleasedSlider( ) +{ + RequestStatisticsUpdate(); +} + +void QmitkPartialVolumeAnalysisView::UpsamplingReleasedSlider( ) +{ + RequestStatisticsUpdate(); +} + +void QmitkPartialVolumeAnalysisView::GaussianSigmaReleasedSlider( ) +{ + RequestStatisticsUpdate(); +} + +void QmitkPartialVolumeAnalysisView::SimilarAnglesReleasedSlider( ) +{ + +} + +void QmitkPartialVolumeAnalysisView::ToClipBoard() +{ + + std::vector* > vals = m_Controls->m_HistogramWidget->m_Vals; + QString clipboardText; + for (std::vector* >::iterator it = vals.begin(); it + != vals.end(); ++it) + { + for (std::vector::iterator it2 = (**it).begin(); it2 != (**it).end(); ++it2) + { + clipboardText.append(QString("%1 \t").arg(*it2)); + } + clipboardText.append(QString("\n")); + } + + QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); +} + +void QmitkPartialVolumeAnalysisView::PropertyChanged(const mitk::DataNode* /*node*/, const mitk::BaseProperty* /*prop*/) +{ +} + +void QmitkPartialVolumeAnalysisView::NodeChanged(const mitk::DataNode* /*node*/) +{ +} + +void QmitkPartialVolumeAnalysisView::NodeRemoved(const mitk::DataNode* node) +{ + + if( node == m_SelectedPlanarFigureNodes->GetNode().GetPointer() + || node == m_SelectedMaskNode.GetPointer() ) + { + this->Select(NULL,true,false); + SetMeasurementInfoToRenderWindow(""); + } + + if( node == m_SelectedImageNodes->GetNode().GetPointer() ) + { + this->Select(NULL,false,true); + SetMeasurementInfoToRenderWindow(""); + } +} + +void QmitkPartialVolumeAnalysisView::NodeAddedInDataStorage(const mitk::DataNode* node) +{ + if(!m_Visible) + return; + + mitk::DataNode* nonConstNode = const_cast(node); + mitk::PlanarFigure* figure = dynamic_cast(nonConstNode->GetData()); + + if(figure) + { + + // set interactor for new node (if not already set) + mitk::PlanarFigureInteractor::Pointer figureInteractor + = dynamic_cast(node->GetInteractor()); + + if(figureInteractor.IsNull()) + figureInteractor = mitk::PlanarFigureInteractor::New("PlanarFigureInteractor", nonConstNode); + + mitk::GlobalInteraction::GetInstance()->AddInteractor(figureInteractor); + + // remove uninitialized old planars + if( m_SelectedPlanarFigureNodes->GetNode().IsNotNull() && m_CurrentFigureNodeInitialized == false ) + { + mitk::Interactor::Pointer oldInteractor = m_SelectedPlanarFigureNodes->GetNode()->GetInteractor(); + if(oldInteractor.IsNotNull()) + mitk::GlobalInteraction::GetInstance()->RemoveInteractor(oldInteractor); + + this->GetDefaultDataStorage()->Remove(m_SelectedPlanarFigureNodes->GetNode()); + } + + } +} + +void QmitkPartialVolumeAnalysisView::TextIntON() +{ + if(m_ClusteringResult.IsNotNull()) + { + if(m_TexIsOn) + { + m_Controls->m_TextureIntON->setIcon(*m_IconTexOFF); + } + else + { + m_Controls->m_TextureIntON->setIcon(*m_IconTexON); + } + + m_ClusteringResult->SetBoolProperty("texture interpolation", !m_TexIsOn); + m_TexIsOn = !m_TexIsOn; + GetActiveStdMultiWidget()->RequestUpdate(); + } +} diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisView.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisView.h new file mode 100644 index 0000000000..f417cb7146 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisView.h @@ -0,0 +1,265 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-05-28 20:08:26 +0200 (Do, 28 Mai 2009) $ +Version: $Revision: 10185 $ + +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. + +=========================================================================*/ + +#if !defined(QmitkPartialVolumeAnalysisView_H__INCLUDED) +#define QmitkPartialVolumeAnalysisView_H__INCLUDED + +#include "QmitkFunctionality.h" +#include "ui_QmitkPartialVolumeAnalysisViewControls.h" + + +// berry +#include +#include + +// itk +#include +#include "itkDiffusionTensorPrincipleDirectionImageFilter.h" + +// qmitk +#include "QmitkStepperAdapter.h" +#include "QmitkRenderWindow.h" + +// mitk +#include "mitkPartialVolumeAnalysisHistogramCalculator.h" +#include "mitkPlanarLine.h" +#include +#include "mitkDataStorageSelection.h" +#include + +// vtk +#include +#include +#include +//#include "itkProcessObject.h" + +/*! +\brief QmitkPartialVolumeAnalysis + +\sa QmitkFunctionality +\ingroup Functionalities +*/ +class QmitkPartialVolumeAnalysisView : public QmitkFunctionality//, public itk::ProcessObject +{ + Q_OBJECT + +public: + + /*! + \ Convenient typedefs + */ + typedef mitk::DataStorage::SetOfObjects ConstVector; + typedef ConstVector::ConstPointer ConstVectorPointer; + typedef ConstVector::ConstIterator ConstVectorIterator; + typedef mitk::PartialVolumeAnalysisHistogramCalculator HistogramCalculatorType; + typedef HistogramCalculatorType::HistogramType HistogramType; + typedef mitk::PartialVolumeAnalysisClusteringCalculator ClusteringType; + typedef itk::DiffusionTensorPrincipleDirectionImageFilter DirectionsFilterType; + + /*! + \brief default constructor + */ + QmitkPartialVolumeAnalysisView(QObject *parent=0, const char *name=0); + + /*! + \brief default destructor + */ + virtual ~QmitkPartialVolumeAnalysisView(); + + /*! + \brief method for creating the widget containing the application controls, like sliders, buttons etc. + */ + virtual void CreateQtPartControl(QWidget *parent); + + /*! + \brief method for creating the connections of main and control widget + */ + virtual void CreateConnections(); + + bool IsExclusiveFunctionality() const; + + virtual bool event( QEvent *event ); + + void OnSelectionChanged( std::vector nodes ); + + virtual void Activated(); + + virtual void Deactivated(); + + bool AssertDrawingIsPossible(bool checked); + + virtual void NodeChanged(const mitk::DataNode* node); + virtual void PropertyChanged(const mitk::DataNode* node, const mitk::BaseProperty* prop); + virtual void NodeRemoved(const mitk::DataNode* node); + virtual void NodeAddedInDataStorage(const mitk::DataNode* node); + + virtual void AddFigureToDataStorage(mitk::PlanarFigure* figure, const QString& name, + const char *propertyKey = NULL, mitk::BaseProperty *property = NULL ); + + void PlanarFigureInitialized(); + + void PlanarFigureFocus(mitk::DataNode* node); + + void ShowClusteringResults(); +protected slots: + + void EstimateCircle(); + void SetHistogramVisibility(); + void SetAdvancedVisibility(); + + void NumberBinsChangedSlider(int v ); + void UpsamplingChangedSlider( int v ); + void GaussianSigmaChangedSlider( int v ); + void SimilarAnglesChangedSlider(int v ); + + void OpacityChangedSlider(int v ); + + void NumberBinsReleasedSlider( ); + void UpsamplingReleasedSlider( ); + void GaussianSigmaReleasedSlider( ); + void SimilarAnglesReleasedSlider( ); + + void ActionDrawEllipseTriggered(); + void ActionDrawRectangleTriggered(); + void ActionDrawPolygonTriggered(); + + void ToClipBoard(); + + void GreenRadio(bool checked); + void PartialVolumeRadio(bool checked); + void BlueRadio(bool checked); + void AllRadio(bool checked); + + void OnRenderWindowDelete(QObject * obj); + + void TextIntON(); + +protected: + + void StdMultiWidgetAvailable( QmitkStdMultiWidget& stdMultiWidget ); + + /** \brief Issues a request to update statistics by sending an event to the + * Qt event processing queue. + * + * Statistics update should only be executed after program execution returns + * to the Qt main loop. This mechanism also prevents multiple execution of + * updates where only one is required.*/ + void RequestStatisticsUpdate(); + + /** \brief Recalculate statistics for currently selected image and mask and + * update the GUI. */ + void UpdateStatistics(); + + /** \brief Listener for progress events to update progress bar. */ + void UpdateProgressBar(); + + /** \brief Removes any cached images which are no longer referenced elsewhere. */ + void RemoveOrphanImages(); + + void Select( mitk::DataNode::Pointer node, bool clearMaskOnFirstArgNULL=false, bool clearImageOnFirstArgNULL=false ); + + void Visible( ); + + void SetMeasurementInfoToRenderWindow(const QString& text); + + void FindRenderWindow(mitk::DataNode* node); + + void ExtractTensorImages( + mitk::Image::ConstPointer tensorimage); + + typedef std::map< mitk::Image *, mitk::PartialVolumeAnalysisHistogramCalculator::Pointer > + PartialVolumeAnalysisMapType; + + /*! + * controls containing sliders for scrolling through the slices + */ + Ui::QmitkPartialVolumeAnalysisViewControls *m_Controls; + + QmitkStepperAdapter* m_TimeStepperAdapter; + unsigned int m_CurrentTime; + + QString m_Clipboard; + + // result text rendering + vtkRenderer * m_MeasurementInfoRenderer; + vtkCornerAnnotation *m_MeasurementInfoAnnotation; + + // Image and mask data + mitk::DataStorageSelection::Pointer m_SelectedImageNodes; + mitk::Image::Pointer m_SelectedImage; + + mitk::DataNode::Pointer m_SelectedMaskNode; + mitk::Image::Pointer m_SelectedImageMask; + + mitk::DataStorageSelection::Pointer m_SelectedPlanarFigureNodes; + mitk::PlanarFigure::Pointer m_SelectedPlanarFigure; + + bool m_IsTensorImage; + mitk::Image::Pointer m_FAImage; + mitk::Image::Pointer m_CAImage; + mitk::Image::Pointer m_RDImage; + mitk::Image::Pointer m_ADImage; + mitk::Image::Pointer m_MDImage; +// mitk::Image::Pointer m_DirectionImage; + mitk::Image::Pointer m_DirectionComp1Image; + mitk::Image::Pointer m_DirectionComp2Image; + mitk::Image::Pointer m_AngularErrorImage; + + QmitkRenderWindow* m_SelectedRenderWindow; + QmitkRenderWindow* m_LastRenderWindow; + + long m_ImageObserverTag; + long m_ImageMaskObserverTag; + long m_PlanarFigureObserverTag; + + // Hash map for associating one image statistics calculator with each iamge + // (so that previously calculated histograms / statistics can be recovered + // if a recalculation is not required) + PartialVolumeAnalysisMapType m_PartialVolumeAnalysisMap; + + HistogramCalculatorType::Pointer m_CurrentStatisticsCalculator; + + bool m_CurrentStatisticsValid; + + bool m_StatisticsUpdatePending; + + bool m_GaussianSigmaChangedSliding; + bool m_NumberBinsSliding; + bool m_UpsamplingChangedSliding; + + mitk::DataNode::Pointer m_ClusteringResult; + + int m_EllipseCounter; + int m_RectangleCounter; + int m_PolygonCounter; + unsigned int m_InitializedObserverTag; + bool m_CurrentFigureNodeInitialized; + + int m_QuantifyClass; + + ClusteringType::HelperStructPerformRGBClusteringRetval* m_CurrentRGBClusteringResults; + ClusteringType::HelperStructPerformClusteringRetval *m_CurrentPerformClusteringResults; + +// mitk::DataNode::Pointer m_newnode; +// mitk::DataNode::Pointer m_newnode2; + QIcon* m_IconTexOFF; + QIcon* m_IconTexON; + bool m_TexIsOn; +}; + + +#endif // !defined(QmitkPartialVolumeAnalysis_H__INCLUDED) diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisViewControls.ui new file mode 100644 index 0000000000..62f1529a98 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPartialVolumeAnalysisViewControls.ui @@ -0,0 +1,762 @@ + + + QmitkPartialVolumeAnalysisViewControls + + + true + + + + 0 + 0 + 327 + 427 + + + + Form + + + + 0 + + + 0 + + + + + false + + + true + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 6 + + + 0 + + + 0 + + + 0 + + + 9 + + + 0 + + + + + + 0 + 0 + + + + Image: + + + + + + + false + + + true + + + + + + + + 0 + 0 + + + + Mask: + + + + + + + false + + + true + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 30 + 30 + + + + + + + + :/QmitkPartialVolumeAnalysisView/circle.png:/QmitkPartialVolumeAnalysisView/circle.png + + + + 32 + 32 + + + + true + + + true + + + + + + + + 30 + 30 + + + + + + + + :/QmitkPartialVolumeAnalysisView/rectangle.png:/QmitkPartialVolumeAnalysisView/rectangle.png + + + + 32 + 32 + + + + true + + + true + + + + + + + + 30 + 30 + + + + + + + + :/QmitkPartialVolumeAnalysisView/polygon.png:/QmitkPartialVolumeAnalysisView/polygon.png + + + + 32 + 32 + + + + true + + + true + + + + + + + + + + + + + true + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + Upsampling + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + 40 + + + 1 + + + 25 + + + Qt::Horizontal + + + + + + + + 50 + 0 + + + + 2.5 + + + + + + + + + + Similar angles + + + + + + + QFrame::NoFrame + + + + 0 + + + + + 90 + + + 0 + + + Qt::Horizontal + + + QSlider::NoTicks + + + + + + + + 50 + 0 + + + + 90° + + + + + + + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + 0 + + + + + display histogram + + + + + + + true + + + + 0 + 0 + + + + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + + 20 + 20 + + + + + + + true + + + true + + + + + + + Green + + + + + + + Partial Volume + + + Partial Volume + + + Partial Volume + + + Partial Volume + + + PV + + + + + + + Red + + + true + + + + + + + All + + + + + + + + + + QFrame::NoFrame + + + + 0 + + + + + Opacity + + + + + + + 10 + + + 5 + + + Qt::Horizontal + + + QSlider::TicksBelow + + + + + + + + + + + + + + + Histogram to Clipboard + + + + + + + + + Advanced + + + + + + + Qt::Vertical + + + QSizePolicy::Preferred + + + + 10 + 1 + + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + QFormLayout::AllNonFixedFieldsGrow + + + 0 + + + + + Blurring + + + + + + + QFrame::NoFrame + + + QFrame::Raised + + + + 0 + + + + + 200 + + + 1 + + + 0 + + + Qt::Horizontal + + + + + + + + 50 + 0 + + + + 0.0 + + + + + + + + + + # Bins + + + + + + + QFrame::NoFrame + + + + 0 + + + + + 100 + + + 10 + + + Qt::Horizontal + + + QSlider::NoTicks + + + + + + + + 50 + 0 + + + + 50 + + + + + + + + + + quantiles + + + + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + 0 + + + + + 1.000000000000000 + + + 0.010000000000000 + + + 0.250000000000000 + + + + + + + 1.000000000000000 + + + 0.010000000000000 + + + 0.750000000000000 + + + + + + + + + + Estimate circle from binary image + + + + + + + + + + Qt::Vertical + + + + 20 + 40 + + + + + + + + + QmitkPartialVolumeAnalysisWidget + QWidget +
QmitkPartialVolumeAnalysisWidget.h
+ 1 + + + + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp index 9d88bb7da8..a95a91f56b 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/mitkPluginActivator.cpp @@ -1,44 +1,46 @@ #include "mitkPluginActivator.h" #include #include "src/internal/QmitkDiffusionImagingPublicPerspective.h" #include "src/internal/QmitkQBallReconstructionView.h" #include "src/internal/QmitkPreprocessingView.h" #include "src/internal/QmitkDiffusionDicomImportView.h" #include "src/internal/QmitkDiffusionQuantificationView.h" #include "src/internal/QmitkTensorReconstructionView.h" #include "src/internal/QmitkControlVisualizationPropertiesView.h" #include "src/internal/QmitkODFDetailsView.h" #include "src/internal/QmitkGlobalFiberTrackingView.h" #include "src/internal/QmitkFiberBundleOperationsView.h" #include "src/internal/QmitkFiberBundleDeveloperView.h" +#include "src/internal/QmitkPartialVolumeAnalysisView.h" namespace mitk { void PluginActivator::start(ctkPluginContext* context) { BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionImagingPublicPerspective, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkQBallReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkPreprocessingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionDicomImport, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkDiffusionQuantificationView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkTensorReconstructionView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkControlVisualizationPropertiesView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkODFDetailsView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkGlobalFiberTrackingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkFiberBundleOperationsView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkFiberBundleDeveloperView, context) + BERRY_REGISTER_EXTENSION_CLASS(QmitkPartialVolumeAnalysisView, context) } void PluginActivator::stop(ctkPluginContext* context) { Q_UNUSED(context) } } Q_EXPORT_PLUGIN2(org_mitk_gui_qt_diffusionimaging, mitk::PluginActivator) diff --git a/Modules/DiffusionImaging/Algorithms/itkCartesianToPolarVectorImageFilter.h b/Modules/DiffusionImaging/Algorithms/itkCartesianToPolarVectorImageFilter.h new file mode 100644 index 0000000000..e298e18d1b --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkCartesianToPolarVectorImageFilter.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 __itkCartesianToPolarVectorImageFilter_h +#define __itkCartesianToPolarVectorImageFilter_h + +#include "itkUnaryFunctorImageFilter.h" + +#define C2P_PI 3.14159265358979323846 + +namespace itk +{ + + namespace Functor { + + template< typename TInput, typename TOutput, bool symmetric > + class CartesianToPolarFunction + { + public: + CartesianToPolarFunction() {} + ~CartesianToPolarFunction() {} + bool operator!=( const CartesianToPolarFunction & ) const + { + return false; + } + bool operator==( const CartesianToPolarFunction & other ) const + { + return !(*this != other); + } + inline TOutput operator()( const TInput & x ) + { + + TOutput opoint; + + if(x[0] || x[1] || x[2]) + { + opoint[0] = sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] ); + opoint[1] = atan2( x[1], x[0] ); + opoint[2] = 0.5*C2P_PI - atan( x[2] / sqrt( x[0] * x[0] + x[1] * x[1] ) ); + + if(symmetric && opoint[1]>C2P_PI) + { + opoint[1] = opoint[1] - C2P_PI; + } + } + else + { + opoint[0] = 0; + opoint[1] = 0; + opoint[2] = 0; + } + return opoint; + + } + }; + + } // end namespace functor + + + /** \class CartesianToPolarVectorImageFilter + * + */ + template + class CartesianToPolarVectorImageFilter : + public + UnaryFunctorImageFilter > + { + public: + /** Standard class typedefs. */ + typedef CartesianToPolarVectorImageFilter Self; + typedef UnaryFunctorImageFilter< + TInputImage,TOutputImage, + Functor::CartesianToPolarFunction< + typename TInputImage::PixelType, + typename TOutputImage::PixelType, symmetric > > Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + typedef typename Superclass::OutputImageType OutputImageType; + typedef typename OutputImageType::PixelType OutputPixelType; + typedef typename TInputImage::PixelType InputPixelType; + typedef typename InputPixelType::ValueType InputValueType; + + /** Run-time type information (and related methods). */ + itkTypeMacro( CartesianToPolarVectorImageFilter, UnaryFunctorImageFilter ); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Print internal ivars */ + void PrintSelf(std::ostream& os, Indent indent) const + { this->Superclass::PrintSelf( os, indent ); } + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(InputHasNumericTraitsCheck, + (Concept::HasNumericTraits)); + /** End concept checking */ +#endif + + protected: + CartesianToPolarVectorImageFilter() {}; + virtual ~CartesianToPolarVectorImageFilter() {}; + + private: + CartesianToPolarVectorImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; + + + +} // end namespace itk + +#endif // __itkCartesianToPolarVectorImageFilter_h diff --git a/Modules/DiffusionImaging/Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h b/Modules/DiffusionImaging/Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h new file mode 100644 index 0000000000..3a6ae37ef7 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h @@ -0,0 +1,89 @@ +/*========================================================================= + + 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 __itkDiffusionTensorPrincipleDirectionImageFilter_h_ +#define __itkDiffusionTensorPrincipleDirectionImageFilter_h_ + +#include "MitkDiffusionImagingMBIExports.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 "itkDiffusionTensor3D.h" + +namespace itk{ + /** \class DiffusionTensorPrincipleDirectionImageFilter + */ + + template< class TTensorPixelType, class TPDPixelType=double> + class DiffusionTensorPrincipleDirectionImageFilter : + public ImageToImageFilter< Image< DiffusionTensor3D, 3 >, + Image< Vector< TPDPixelType, 3 >, 3 > > + { + + public: + + typedef DiffusionTensorPrincipleDirectionImageFilter Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ImageToImageFilter< Image< DiffusionTensor3D, 3 >, + Image< Vector< TPDPixelType, 3 >, 3 > > + Superclass; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(DiffusionTensorPrincipleDirectionImageFilter, + ImageToImageFilter); + + typedef TTensorPixelType TensorComponentType; + + typedef TPDPixelType DirectionPixelType; + + typedef typename Superclass::InputImageType InputImageType; + + typedef typename Superclass::OutputImageType OutputImageType; + + typedef typename Superclass::OutputImageRegionType + OutputImageRegionType; + + void SetImage( const InputImageType *image ); + + protected: + DiffusionTensorPrincipleDirectionImageFilter(); + ~DiffusionTensorPrincipleDirectionImageFilter() {}; + void PrintSelf(std::ostream& os, Indent indent) const; + + void BeforeThreadedGenerateData(); + void ThreadedGenerateData( const + OutputImageRegionType &outputRegionForThread, int); + + private: + + }; + +} + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkDiffusionTensorPrincipleDirectionImageFilter.txx" +#endif + +#endif //__itkDiffusionTensorPrincipleDirectionImageFilter_h_ + diff --git a/Modules/DiffusionImaging/Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.txx b/Modules/DiffusionImaging/Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.txx new file mode 100644 index 0000000000..f759c45479 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.txx @@ -0,0 +1,118 @@ + +#ifndef __itkDiffusionTensorPrincipleDirectionImageFilter_txx +#define __itkDiffusionTensorPrincipleDirectionImageFilter_txx + +#include +#include +#include + +//#include "itkDiffusionTensorPrincipleDirectionImageFilter.h" +#include "itkImageRegionConstIterator.h" +#include "itkImageRegionConstIteratorWithIndex.h" +#include "itkImageRegionIterator.h" +#include "itkArray.h" +#include "vnl/vnl_vector.h" + +namespace itk { + + //#define QBALL_RECON_PI 3.14159265358979323846 + + template< class TTensorPixelType, class TPDPixelType> + DiffusionTensorPrincipleDirectionImageFilter< TTensorPixelType, + TPDPixelType> + ::DiffusionTensorPrincipleDirectionImageFilter() + { + // At least 1 inputs is necessary for a vector image. + // For images added one at a time we need at least six + this->SetNumberOfRequiredInputs( 1 ); + } + + template< class TTensorPixelType, + class TPDPixelType> + void DiffusionTensorPrincipleDirectionImageFilter< TTensorPixelType, + TPDPixelType> + ::BeforeThreadedGenerateData() + { + } + + template< class TTensorPixelType, + class TPDPixelType> + void DiffusionTensorPrincipleDirectionImageFilter< TTensorPixelType, + TPDPixelType> + ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, + int ) + { + + typedef itk::DiffusionTensor3D TensorType; + typedef ImageRegionConstIterator< InputImageType > InputIteratorType; + typedef typename InputImageType::PixelType InputTensorType; + typename InputImageType::Pointer inputImagePointer = NULL; + inputImagePointer = static_cast< InputImageType * >( + this->ProcessObject::GetInput(0) ); + + typename OutputImageType::Pointer outputImage = + static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); + outputImage->SetSpacing( inputImagePointer->GetSpacing() ); // Set the image spacing + outputImage->SetOrigin( inputImagePointer->GetOrigin() ); // Set the image origin + outputImage->SetDirection( inputImagePointer->GetDirection() ); // Set the image direction + outputImage->SetRegions( inputImagePointer->GetLargestPossibleRegion() ); + outputImage->Allocate(); + + ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); + oit.GoToBegin(); + + InputIteratorType git(inputImagePointer, outputRegionForThread ); + git.GoToBegin(); + while( !git.IsAtEnd() ) + { + InputTensorType b = git.Get(); + TensorType tensor = b.GetDataPointer(); + + typename OutputImageType::PixelType dir; + typename TensorType::EigenValuesArrayType eigenvalues; + typename TensorType::EigenVectorsMatrixType eigenvectors; + if(tensor.GetTrace()!=0) + { + tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); + +// int index = 2; +// if( (eigenvalues[0] >= eigenvalues[1]) +// && (eigenvalues[0] >= eigenvalues[2]) ) +// index = 0; +// else if(eigenvalues[1] >= eigenvalues[2]) +// index = 1; + + vnl_vector_fixed vec; + vec[0] = eigenvectors(/*index*/2,0); + vec[1] = eigenvectors(/*index*/2,1); + vec[2] = eigenvectors(/*index*/2,2); + + dir[0] = (TPDPixelType)vec[0]; + dir[1] = (TPDPixelType)vec[1]; + dir[2] = (TPDPixelType)vec[2]; + } + else + { + dir[0] = (TPDPixelType)0; + dir[1] = (TPDPixelType)0; + dir[2] = (TPDPixelType)0; + } + oit.Set( dir ); + + ++oit; + ++git; // Gradient image iterator + } + + std::cout << "One Thread finished extraction" << std::endl; + } + + template< class TTensorPixelType, + class TPDPixelType> + void DiffusionTensorPrincipleDirectionImageFilter< TTensorPixelType, + TPDPixelType> + ::PrintSelf(std::ostream& os, Indent indent) const + { + } + +} +#endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/Algorithms/itkPolarToCartesianVectorImageFilter.h b/Modules/DiffusionImaging/Algorithms/itkPolarToCartesianVectorImageFilter.h new file mode 100644 index 0000000000..a9ccb5eabe --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkPolarToCartesianVectorImageFilter.h @@ -0,0 +1,123 @@ +/*========================================================================= + +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 __itkPolarToCartesianVectorImageFilter_h +#define __itkPolarToCartesianVectorImageFilter_h + +#include "itkUnaryFunctorImageFilter.h" + +#define P2C_PI 3.14159265358979323846 + +namespace itk +{ + + namespace Functor { + + template< typename TInput, typename TOutput, bool symmetric > + class PolarToCartesianFunction + { + public: + PolarToCartesianFunction() {} + ~PolarToCartesianFunction() {} + bool operator!=( const PolarToCartesianFunction & ) const + { + return false; + } + bool operator==( const PolarToCartesianFunction & other ) const + { + return !(*this != other); + } + inline TOutput operator()( const TInput & x ) + { + + TOutput opoint; + opoint[0] = x[0] * cos( x[1] ) * sin( x[2] ); + opoint[1] = x[0] * sin( x[1] ) * sin( x[2] ); + opoint[2] = x[0] * cos( x[2] ); + + if(symmetric && opoint[2]<0) + { + opoint[2] = -opoint[2]; + } + + return opoint; + ; + } + }; + + } // end namespace functor + + + /** \class PolarToCartesianVectorImageFilter + * + */ + template + class PolarToCartesianVectorImageFilter : + public + UnaryFunctorImageFilter > + { + public: + /** Standard class typedefs. */ + typedef PolarToCartesianVectorImageFilter Self; + typedef UnaryFunctorImageFilter< + TInputImage,TOutputImage, + Functor::PolarToCartesianFunction< + typename TInputImage::PixelType, + typename TOutputImage::PixelType, symmetric > > Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + typedef typename Superclass::OutputImageType OutputImageType; + typedef typename OutputImageType::PixelType OutputPixelType; + typedef typename TInputImage::PixelType InputPixelType; + typedef typename InputPixelType::ValueType InputValueType; + + /** Run-time type information (and related methods). */ + itkTypeMacro( PolarToCartesianVectorImageFilter, UnaryFunctorImageFilter ); + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Print internal ivars */ + void PrintSelf(std::ostream& os, Indent indent) const + { this->Superclass::PrintSelf( os, indent ); } + +#ifdef ITK_USE_CONCEPT_CHECKING + /** Begin concept checking */ + itkConceptMacro(InputHasNumericTraitsCheck, + (Concept::HasNumericTraits)); + /** End concept checking */ +#endif + + protected: + PolarToCartesianVectorImageFilter() {}; + virtual ~PolarToCartesianVectorImageFilter() {}; + + private: + PolarToCartesianVectorImageFilter(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; + + + +} // end namespace itk + +#endif // __itkPolarToCartesianVectorImageFilter_h diff --git a/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp new file mode 100644 index 0000000000..56bed42c5f --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp @@ -0,0 +1,987 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2009-05-12 19:56:03 +0200 (Di, 12 Mai 2009) $ +Version: $Revision: 17179 $ + +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 "mitkPartialVolumeAnalysisClusteringCalculator.h" + +#include "mitkImageAccessByItk.h" +#include "mitkImageToItk.h" +#include "itkScalarImageToHistogramGenerator.h" +#include "itkListSample.h" + +#define PVA_PI 3.14159265358979323846 + +namespace mitk +{ + + PartialVolumeAnalysisClusteringCalculator::PartialVolumeAnalysisClusteringCalculator() + : m_MaxIt(100), m_StepsNumIntegration(100) + { + } + + + PartialVolumeAnalysisClusteringCalculator::~PartialVolumeAnalysisClusteringCalculator() + { + } + + PartialVolumeAnalysisClusteringCalculator::ParamsType* + PartialVolumeAnalysisClusteringCalculator::InitialGuess(HistType h) const + { + int s = h.xVals.size(); + int s30 = 0.3 * s; + int s70 = 0.7 * s; + ParamsType* params = new ParamsType(); + params->means[0] = h.xVals(s30); + params->means[1] = h.xVals(s70); + params->sigmas[0] = (h.xVals(s-1) - h.xVals(0)) / 5.0; + params->sigmas[1] = (h.xVals(s-1) - h.xVals(0)) / 5.0; + params->ps[0] = 0.4; + params->ps[1] = 0.4; + + return params; + } + + PartialVolumeAnalysisClusteringCalculator::ParamsType* + PartialVolumeAnalysisClusteringCalculator::Cluster(HistType h, ParamsType *initialGuess) const + { + ParamsType *params = new ParamsType(); + params->Initialize(initialGuess); + + double ll = 9999999999999999999.9; + double oll; + double sml = (h.xVals[1] - h.xVals[0]) / 1000; + + int arraysz = h.xVals.size(); + + for (unsigned int it = 0; it < m_MaxIt; it++) + { + // wie sehen basisfunktionen zu aktuellen params aus? + ClusterResultType curves = CalculateCurves(*params,h.xVals); + + // summe der basisfunktionen + for(int i=0; i2 && oll-ll < arraysz/2*1e-15 ) + { + break; + } + + for(int j=0; j<2; j++) + { + VecType array, p(arraysz); + array = curves.vals[j]; + + for(int i=0; ips[j] = 0; + params->means[j] = 0; + for(int i=0; ips[j] += p(i); + params->means[j] += h.xVals(i)*p(i); + } + params->means[j] /= params->ps[j]; + + VecType vr = h.xVals; + for(int i=0; imeans[j]; + } + + params->sigmas[j] = 0; + for(int i=0; isigmas[j] += vr(i)*vr(i)*p(i); + } + params->sigmas[j] /= params->ps[j]; + params->sigmas[j] += sml; + + } + + double p3 = 0; + for(int i=0; ips[j] = params->ps[j] + 1e-3; + sum += params->ps[j]; + } + sum += p3; + + for(int j=0; j<2; j++) + { + params->ps[j] = params->ps[j] / sum; + } + p3 /= sum; + + } + + return params; + + } + + void PartialVolumeAnalysisClusteringCalculator::Normalize(ParamsType params, ClusterResultType* curves) const + { + double sum1=0, sum2=0, sum3=0; + int arraysz = curves->vals[0].size(); + + for(int i=0; ivals[0](i); + sum2 += curves->vals[1](i); + sum3 += curves->mixedVals[0](i); + } + + sum1 /= params.ps[0]; + sum2 /= params.ps[1]; + sum3 /= 1.0 -params.ps[0] -params.ps[1]; + + for(int i=0; ivals[0](i) /= sum1; + curves->vals[1](i) /= sum2; + curves->mixedVals[0](i) /= sum3; + } + + for(int i=0; icombiVals(i) = curves->mixedVals[0](i) + curves->vals[0](i) + curves->vals[1](i); + } + } + + PartialVolumeAnalysisClusteringCalculator::ClusterResultType + PartialVolumeAnalysisClusteringCalculator::CalculateCurves(ParamsType params, VecType xVals) const + { + + int arraysz = xVals.size(); + ClusterResultType result(arraysz); + + for( int j=0; j<2; j++) + { + for(int i=0; i 0.0000000000000001) + { + itprob.Set( aposteriori ); + maxp = aposteriori > maxp ? aposteriori : maxp; + } + else + { + itprob.Set(0.0f); + } + } + + ++itimage; + ++itprob; + } + + itprob = itprob.Begin(); + itdisp = itdisp.Begin(); + + while( !itprob.IsAtEnd() ) + { + if(itprob.Get()) + { + typename DisplayImageType::PixelType rgba; + rgba.Set(255.0f, 0.0f, 0.0f, 255.0f*(itprob.Get()/maxp)); + itdisp.Set( rgba ); + } + ++itprob; + ++itdisp; + } + + outImage1->InitializeByItk(probimage.GetPointer()); + outImage1->SetVolume(probimage->GetBufferPointer()); + + outImage2->InitializeByItk(displayimage.GetPointer()); + outImage2->SetVolume(displayimage->GetBufferPointer()); + + } + + + double* PartialVolumeAnalysisClusteringCalculator::PerformQuantification( + mitk::Image::ConstPointer image, mitk::Image::Pointer clusteredImage, mitk::Image::Pointer mask) const + { + + double *retval = new double[2]; + + AccessFixedDimensionByItk_3( + image.GetPointer(), + InternalQuantify, + 3, + clusteredImage.GetPointer(), + retval, mask ); + + return retval; + + } + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisClusteringCalculator::InternalQuantify( + const itk::Image< TPixel, VImageDimension > *image, + mitk::Image::ConstPointer clusteredImage, double* retval, mitk::Image::Pointer mask ) const + { + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< float, VImageDimension > ProbImageType; + typedef itk::Image< unsigned char, VImageDimension > MaskImageType; + + typedef mitk::ImageToItk CastFilterType; + typename CastFilterType::Pointer castFilter = CastFilterType::New(); + castFilter->SetInput( clusteredImage ); + castFilter->Update(); + typename ProbImageType::Pointer clusterImage = castFilter->GetOutput(); + + typename MaskImageType::Pointer itkmask = 0; + if(mask.IsNotNull()) + { + typedef mitk::ImageToItk CastFilterType2; + typename CastFilterType2::Pointer castFilter2 = CastFilterType2::New(); + castFilter2->SetInput( mask ); + castFilter2->Update(); + itkmask = castFilter2->GetOutput(); + } + else + { + itkmask = MaskImageType::New(); + itkmask->SetSpacing( clusterImage->GetSpacing() ); // Set the image spacing + itkmask->SetOrigin( clusterImage->GetOrigin() ); // Set the image origin + itkmask->SetDirection( clusterImage->GetDirection() ); // Set the image direction + itkmask->SetRegions( clusterImage->GetLargestPossibleRegion() ); + itkmask->Allocate(); + itkmask->FillBuffer(1); + } + + itk::ImageRegionConstIterator + itimage(image, image->GetLargestPossibleRegion()); + + itk::ImageRegionConstIterator + itprob(clusterImage, clusterImage->GetLargestPossibleRegion()); + + itk::ImageRegionConstIterator + itmask(itkmask, itkmask->GetLargestPossibleRegion()); + + itimage = itimage.Begin(); + itprob = itprob.Begin(); + itmask = itmask.Begin(); + + double totalProb = 0; + double measurement = 0; + double error = 0; + + while( !itimage.IsAtEnd() && !itprob.IsAtEnd() && !itmask.IsAtEnd() ) + { + double valImag = itimage.Get(); + double valProb = itprob.Get(); + double valMask = itmask.Get(); + + typename ProbImageType::PixelType prop = valProb * valMask; + + totalProb += prop; + measurement += valImag * prop; + error += valImag * valImag * prop; + + ++itimage; + ++itprob; + ++itmask; + } + + measurement = measurement / totalProb; + error = error / totalProb; + retval[0] = measurement; + retval[1] = sqrt( error - measurement*measurement ); + + } + + + mitk::Image::Pointer PartialVolumeAnalysisClusteringCalculator::CaculateAngularErrorImage( + mitk::Image::Pointer comp1, mitk::Image::Pointer comp2, mitk::Image::Pointer probImg) const + { + + // cast input images to itk + typedef itk::Image ImageType; + typedef mitk::ImageToItk CastType; + CastType::Pointer caster = CastType::New(); + caster->SetInput(comp1); + caster->Update(); + ImageType::Pointer comp1Image = caster->GetOutput(); + + caster = CastType::New(); + caster->SetInput(comp2); + caster->Update(); + ImageType::Pointer comp2Image = caster->GetOutput(); + + caster = CastType::New(); + caster->SetInput(probImg); + caster->Update(); + ImageType::Pointer probImage = caster->GetOutput(); + + // figure out maximum probability for fiber class + float maxProb = 0; + itk::ImageRegionConstIterator + itprob(probImage, probImage->GetLargestPossibleRegion()); + itprob = itprob.Begin(); + while( !itprob.IsAtEnd() ) + { + maxProb = itprob.Get() > maxProb ? itprob.Get() : maxProb; + ++itprob; + } + + // generate a list sample of angles at positions + // where the fiber-prob is higher than .2*maxprob + typedef float MeasurementType; + const unsigned int MeasurementVectorLength = 2; + typedef itk::Vector< MeasurementType , MeasurementVectorLength > + MeasurementVectorType; + typedef itk::Statistics::ListSample< MeasurementVectorType > ListSampleType; + ListSampleType::Pointer listSample = ListSampleType::New(); + listSample->SetMeasurementVectorSize( MeasurementVectorLength ); + + itk::ImageRegionIterator + it1(comp1Image, comp1Image->GetLargestPossibleRegion()); + itk::ImageRegionIterator + it2(comp2Image, comp2Image->GetLargestPossibleRegion()); + + it1 = it1.Begin(); + it2 = it2.Begin(); + itprob = itprob.Begin(); + + while( !itprob.IsAtEnd() ) + { + if(itprob.Get() > 0.2 * maxProb) + { + MeasurementVectorType mv; + mv[0] = ( MeasurementType ) it1.Get(); + mv[1] = ( MeasurementType ) it2.Get(); + listSample->PushBack(mv); + } + ++it1; + ++it2; + ++itprob; + } + + // generate a histogram from the list sample + typedef float HistogramMeasurementType; + typedef itk::Statistics::ListSampleToHistogramGenerator + < ListSampleType, HistogramMeasurementType, + itk::Statistics::DenseFrequencyContainer, + MeasurementVectorLength > GeneratorType; + GeneratorType::Pointer generator = GeneratorType::New(); + + GeneratorType::HistogramType::SizeType size; + size.Fill(30); + generator->SetNumberOfBins( size ); + + generator->SetListSample( listSample ); + generator->SetMarginalScale( 10.0 ); + + MeasurementVectorType min; + min[0] = ( MeasurementType ) 0; + min[1] = ( MeasurementType ) 0; + generator->SetHistogramMin(min); + + MeasurementVectorType max; + max[0] = ( MeasurementType ) PVA_PI; + max[1] = ( MeasurementType ) PVA_PI; + generator->SetHistogramMax(max); + + generator->Update(); + + // look for frequency mode in the histogram + GeneratorType::HistogramType::ConstPointer histogram = generator->GetOutput(); + GeneratorType::HistogramType::ConstIterator iter = histogram->Begin(); + float maxFreq = 0; + MeasurementVectorType maxValue; + while ( iter != histogram->End() ) + { + if(iter.GetFrequency() > maxFreq) + { + maxFreq = iter.GetFrequency(); + maxValue[0] = iter.GetMeasurementVector()[0]; + maxValue[1] = iter.GetMeasurementVector()[1]; + } + ++iter; + } + + // generate return image that contains the angular + // error of the voxels to the histogram max measurement + ImageType::Pointer returnImage = ImageType::New(); + returnImage->SetSpacing( comp1Image->GetSpacing() ); // Set the image spacing + returnImage->SetOrigin( comp1Image->GetOrigin() ); // Set the image origin + returnImage->SetDirection( comp1Image->GetDirection() ); // Set the image direction + returnImage->SetRegions( comp1Image->GetLargestPossibleRegion() ); + returnImage->Allocate(); + + itk::ImageRegionConstIterator + cit1(comp1Image, comp1Image->GetLargestPossibleRegion()); + + itk::ImageRegionConstIterator + cit2(comp2Image, comp2Image->GetLargestPossibleRegion()); + + itk::ImageRegionIterator + itout(returnImage, returnImage->GetLargestPossibleRegion()); + + cit1 = cit1.Begin(); + cit2 = cit2.Begin(); + itout = itout.Begin(); + + vnl_vector v(3); + v[0] = cos( maxValue[0] ) * sin( maxValue[1] ); + v[1] = sin( maxValue[0] ) * sin( maxValue[1] ); + v[2] = cos( maxValue[1] ); +// MITK_INFO << "max vector: " << v; + while( !cit1.IsAtEnd() ) + { + vnl_vector v1(3); + v1[0] = cos( cit1.Get() ) * sin( cit2.Get() ); + v1[1] = sin( cit1.Get() ) * sin( cit2.Get() ); + v1[2] = cos( cit2.Get() ); + + itout.Set(fabs(angle(v,v1))); +// MITK_INFO << "ang_error " << v1 << ": " << fabs(angle(v,v1)); + + ++cit1; + ++cit2; + ++itout; + } + + mitk::Image::Pointer retval = mitk::Image::New(); + retval->InitializeByItk(returnImage.GetPointer()); + retval->SetVolume(returnImage->GetBufferPointer()); + return retval; + } + + PartialVolumeAnalysisClusteringCalculator::HelperStructPerformRGBClusteringRetval* + PartialVolumeAnalysisClusteringCalculator::PerformRGBQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2) const + { + + HelperStructRGBChannels *rgbChannels = new HelperStructRGBChannels(); + + HelperStructPerformClusteringRetval *resultr = PerformQuantiles(image, histogram, p2, 999999999.0 ); + rgbChannels->r = resultr->clusteredImage; + + HelperStructPerformClusteringRetval *resultg = PerformQuantiles(image, histogram, -999999999.0, p1 ); + rgbChannels->g = resultg->clusteredImage; + + HelperStructPerformClusteringRetval *resultb = PerformQuantiles(image, histogram, p1, p2 ); + rgbChannels->b = resultb->clusteredImage; + + mitk::Image::Pointer outImage = mitk::Image::New(); + + AccessFixedDimensionByItk_2( + rgbChannels->r.GetPointer(), + InternalGenerateRGB, + 3, + rgbChannels, + outImage); + + HelperStructPerformRGBClusteringRetval *retval + = new HelperStructPerformRGBClusteringRetval(); + retval->rgbChannels = rgbChannels; + retval->rgb = outImage; + retval->params = resultr->params; + retval->result = resultr->result; + retval->hist = resultr->hist; + + delete resultr; + delete resultg; + + return retval; + + } + + PartialVolumeAnalysisClusteringCalculator::HelperStructPerformClusteringRetval* + PartialVolumeAnalysisClusteringCalculator::PerformQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2 ) const + { + + HelperStructPerformClusteringRetval *retval = + new HelperStructPerformClusteringRetval(); + + retval->hist = new HistType(); + retval->hist->InitByMitkHistogram(histogram); + + double sum = histogram->GetTotalFrequency(); + + double* q = new double[2]; + q[0] = histogram->Quantile(0, p1); + q[1] = histogram->Quantile(0, p2); + + mitk::Image::Pointer outImage1 = mitk::Image::New(); + mitk::Image::Pointer outImage2 = mitk::Image::New(); + + AccessFixedDimensionByItk_3( + image.GetPointer(), + InternalGenerateQuantileImage, + 3, q, + outImage1, outImage2); + + retval->clusteredImage = outImage1; + retval->displayImage = outImage2; + + delete[] q; + return retval; + + } + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisClusteringCalculator::InternalGenerateQuantileImage( + const itk::Image< TPixel, VImageDimension > *image, double* q, + mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2 ) const + { + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< itk::RGBAPixel, VImageDimension > DisplayImageType; + typedef itk::Image< float, VImageDimension > ProbImageType; + + typename ProbImageType::Pointer probimage = ProbImageType::New(); + probimage->SetSpacing( image->GetSpacing() ); // Set the image spacing + probimage->SetOrigin( image->GetOrigin() ); // Set the image origin + probimage->SetDirection( image->GetDirection() ); // Set the image direction + probimage->SetRegions( image->GetLargestPossibleRegion() ); + probimage->Allocate(); + probimage->FillBuffer(0); + + typename DisplayImageType::Pointer displayimage = DisplayImageType::New(); + displayimage->SetSpacing( image->GetSpacing() ); // Set the image spacing + displayimage->SetOrigin( image->GetOrigin() ); // Set the image origin + displayimage->SetDirection( image->GetDirection() ); // Set the image direction + displayimage->SetRegions( image->GetLargestPossibleRegion() ); + displayimage->Allocate(); + + typename DisplayImageType::PixelType rgba; + rgba.Set(0.0f, 0.0f, 0.0f, 0.0f); + displayimage->FillBuffer(rgba); + + itk::ImageRegionConstIterator + itimage(image, image->GetLargestPossibleRegion()); + + itk::ImageRegionIterator + itprob(probimage, probimage->GetLargestPossibleRegion()); + + itk::ImageRegionIterator + itdisp(displayimage, displayimage->GetLargestPossibleRegion()); + + itimage = itimage.Begin(); + itprob = itprob.Begin(); + + while( !itimage.IsAtEnd() ) + { + if(itimage.Get() > q[0] && itimage.Get() < q[1]) + { + itprob.Set(1.0f); + } + + ++itimage; + ++itprob; + } + + itprob = itprob.Begin(); + itdisp = itdisp.Begin(); + + while( !itprob.IsAtEnd() ) + { + if(itprob.Get()) + { + typename DisplayImageType::PixelType rgba; + rgba.Set(255.0f, 0.0f, 0.0f, 255.0f); + itdisp.Set( rgba ); + } + ++itprob; + ++itdisp; + } + + outImage1->InitializeByItk(probimage.GetPointer()); + outImage1->SetVolume(probimage->GetBufferPointer()); + + outImage2->InitializeByItk(displayimage.GetPointer()); + outImage2->SetVolume(displayimage->GetBufferPointer()); + + } + +} diff --git a/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h new file mode 100644 index 0000000000..ff1040ffc9 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h @@ -0,0 +1,521 @@ +/*======================================================================== +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + + +#ifndef _MITK_PartialVolumeAnalysisClusteringCalculator_H +#define _MITK_PartialVolumeAnalysisClusteringCalculator_H + +#include "MitkDiffusionImagingExports.h" + +#include "mitkCommon.h" +#include "mitkImage.h" + +#include +#include +#include +#include + +#include + +namespace mitk +{ + + class MitkDiffusionImaging_EXPORT PartialVolumeAnalysisClusteringCalculator : public itk::Object + { + public: + + typedef vnl_vector VecType; + typedef mitk::Image::HistogramType MitkHistType; + + class HistType + { + public: + + HistType() + { + } + + HistType(const HistType& p) + { + this->xVals = p.xVals; + this->hVals = p.hVals; + } + + ~HistType() + { + } + + HistType operator= (const HistType& p) + { + if (this != &p) // protect against invalid self-assignment + { + this->xVals = p.xVals; + this->hVals = p.hVals; + } + return *this; + } + + void InitByMitkHistogram(const MitkHistType* histogram) + { + xVals.set_size(histogram->GetSize()[0]); + hVals.set_size(histogram->GetSize()[0]); + double sum = histogram->GetTotalFrequency(); + MitkHistType::ConstIterator endIt = histogram->End(); + MitkHistType::ConstIterator it; + bool firstNonEmptyBinFound = false; + it = histogram->Begin(); + //++it; + int i=0; + while (it != endIt) + { + if(it.GetFrequency() || firstNonEmptyBinFound) + { + firstNonEmptyBinFound = true; + xVals(i) = it.GetMeasurementVector().GetElement(0); + hVals(i) = it.GetFrequency()/sum; + } + ++i; + ++it; + } + + } + + void SetXVals(VecType x) + {xVals = x;} + + void SetHVals(VecType h) + {hVals = h;} + + void SetXVals(std::vector x) + { + int s = x.size(); + xVals.set_size(s); + for(int i=0; i h) + { + int s = h.size(); + hVals.set_size(s); + for(int i=0; i* GetXVals() + { + int s = xVals.size(); + std::vector* retval = new std::vector(s); + for(int i=0; i* GetHVals() + { + int s = hVals.size(); + std::vector* retval = new std::vector(s); + for(int i=0; ivals.clear(); +// this->mixedVals.clear(); + +// int s = p.vals.size(); +// for(int i=0; ivals.push_back(v); +// } + +// s = p.mixedVals.size(); +// for(int i=0; imixedVals.push_back(v); +// } + +// this->combiVals = p.combiVals; +// } +// return *this; +// } + + void Initialize (const ClusterResultType *p) + { + if (this != p) // protect against invalid self-assignment + { + this->vals.clear(); + this->mixedVals.clear(); + + int s = p->vals.size(); + for(int i=0; ivals[i]; + this->vals.push_back(v); + } + + s = p->mixedVals.size(); + for(int i=0; imixedVals[i]; + this->mixedVals.push_back(v); + } + + this->combiVals = p->combiVals; + } + } + + std::vector GetFiberVals() + { + if(vals.size()==2 && vals[0].data_block()) + { + int s = vals[0].size(); + std::vector retval(s); + for(int i=0; i(0); + } + } + + std::vector GetNonFiberVals() + { + if(vals.size()==2 && vals[1].data_block()) + { + int s = vals[1].size(); + std::vector retval(s); + for(int i=0; i(0); + } + } + + std::vector GetMixedVals() + { + if(mixedVals.size()==1 && mixedVals[0].data_block()) + { + int s = mixedVals[0].size(); + std::vector retval(s); + for(int i=0; i(0); + } + } + + std::vector GetCombiVals() + { + if(combiVals.data_block()) + { + int s = combiVals.size(); + std::vector retval(s); + for(int i=0; i(0); + } + } + + + void Print(VecType vec, int nr=10) + { + int sz = vec.size(); + int incr = (int)((1.0*sz)/(1.0*nr)); + for(int i=0; i vals; + std::vector mixedVals; + VecType combiVals; + + }; + + class ParamsType + { + public: + ParamsType() + { + } + + ~ParamsType() + { + } + + void Initialize(const ParamsType *p) + { + if (this != p) // protect against invalid self-assignment + { + means[0] = p->means[0]; + means[1] = p->means[1]; + sigmas[0] = p->sigmas[0]; + sigmas[1] = p->sigmas[1]; + ps[0] = p->ps[0]; + ps[1] = p->ps[1]; + } + } + + void Print() + { + MITK_INFO << "PARAMS" << std::endl; + MITK_INFO << "Class 1: " << means[0] << " +- " << sqrt(sigmas[0]) << " (p=" << ps[0] << ")" << std::endl; + MITK_INFO << "Class 2: " << means[1] << " +- " << sqrt(sigmas[1]) << " (p=" << ps[1] << ")" << std::endl; + MITK_INFO << "Partial V: p=" << 1.0-ps[0]-ps[1] << std::endl; + } + + double means[2]; + double sigmas[2]; + double ps[2]; + }; + + struct HelperStructClusteringResults + { + MitkHistType *interestingHist; + MitkHistType *totalHist; + double p_interesting; + }; + + struct HelperStructRGBChannels + { + mitk::Image::Pointer r; + mitk::Image::Pointer g; + mitk::Image::Pointer b; + + ~HelperStructRGBChannels() + { + r = 0; + g = 0; + b = 0; + } + }; + + struct HelperStructPerformRGBClusteringRetval + { + HelperStructRGBChannels *rgbChannels; + mitk::Image::Pointer rgb; + + ParamsType *params; + ClusterResultType *result; + HistType *hist; + + HelperStructPerformRGBClusteringRetval() : + rgbChannels(0), params(0), result(0), hist(0) + { + } + + ~HelperStructPerformRGBClusteringRetval() + { + rgb = 0; + delete rgbChannels; + } + }; + + struct HelperStructPerformClusteringRetval + { + mitk::Image::Pointer clusteredImage; + mitk::Image::Pointer displayImage; + + ParamsType *params; + ClusterResultType *result; + HistType *hist; + + HelperStructPerformClusteringRetval() : + clusteredImage(0), displayImage(0), + params(0), result(0), hist(0) + { + } + + ~HelperStructPerformClusteringRetval() + { + clusteredImage = 0; + displayImage = 0; + } + }; + + mitkClassMacro( PartialVolumeAnalysisClusteringCalculator, itk::Object ) + itkNewMacro( PartialVolumeAnalysisClusteringCalculator ) + + ParamsType *InitialGuess(HistType h) const; + + ParamsType *Cluster(const HistType h, ParamsType* initialGuess) const; + + ClusterResultType CalculateCurves(ParamsType params, VecType xVals) const; + + void Normalize(ParamsType params, ClusterResultType* curves) const; + + HelperStructPerformClusteringRetval* PerformClustering(mitk::Image::ConstPointer image, const MitkHistType *histogram, int classIdent, HelperStructPerformClusteringRetval* precResult = 0) const; + + HelperStructPerformRGBClusteringRetval* PerformRGBClustering(mitk::Image::ConstPointer image, const MitkHistType *histogram) const; + + HelperStructPerformClusteringRetval* PerformQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2) const; + + HelperStructPerformRGBClusteringRetval* PerformRGBQuantiles(mitk::Image::ConstPointer image, const MitkHistType *histogram, double p1, double p2 ) const; + + double* PerformQuantification(mitk::Image::ConstPointer image, mitk::Image::Pointer clusteredImage, mitk::Image::Pointer mask = 0) const; + + mitk::Image::Pointer CaculateAngularErrorImage( + mitk::Image::Pointer comp1, mitk::Image::Pointer comp2, mitk::Image::Pointer probImg) const; + + template < typename TPixel, unsigned int VImageDimension > + void InternalGenerateRGB( + const itk::Image< TPixel, VImageDimension > *image, + HelperStructRGBChannels *rgb, mitk::Image::Pointer retval ) const; + + template < typename TPixel, unsigned int VImageDimension > + void InternalGenerateProbabilityImage( + const itk::Image< TPixel, VImageDimension > *image, + const HelperStructClusteringResults clusterResults, + mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2 ) const; + + template < typename TPixel, unsigned int VImageDimension > + void InternalQuantify( + const itk::Image< TPixel, VImageDimension > *image, + mitk::Image::ConstPointer clusteredImage, double* retval, mitk::Image::Pointer mask ) const; + + template < typename TPixel, unsigned int VImageDimension > + void InternalGenerateQuantileImage( + const itk::Image< TPixel, VImageDimension > *image, + double* q, + mitk::Image::Pointer outImage1, mitk::Image::Pointer outImage2 ) const; + + ParamsType* Cluster(const HistType h) const + {return Cluster(h, InitialGuess(h));} + + void SetMaxIt(unsigned int it) + { m_MaxIt = it; } + + unsigned int GetMaxIt() + { return m_MaxIt; } + + void SetStepsNumIntegration(unsigned int n) + { m_StepsNumIntegration = n; } + + unsigned int GetStepsNumIntegration() + { return m_StepsNumIntegration; } + + protected: + + PartialVolumeAnalysisClusteringCalculator(); + + virtual ~PartialVolumeAnalysisClusteringCalculator(); + + unsigned int m_MaxIt; + unsigned int m_StepsNumIntegration; + + }; + +} + +#endif // #define _MITK_PartialVolumeAnalysisClusteringCalculator_H diff --git a/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp new file mode 100644 index 0000000000..c08539774d --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp @@ -0,0 +1,1212 @@ +/*========================================================================= +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 "mitkPartialVolumeAnalysisHistogramCalculator.h" +#include "mitkImageAccessByItk.h" +#include "mitkExtractImageFilter.h" + +#include + +#include +#include +#include +#include +#include "itkResampleImageFilter.h" +#include "itkGaussianInterpolateImageFunction.h" + +#include +#include +#include +#include + + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include "itkGaussianInterpolateImageFunction.h" +#include "itkBSplineInterpolateImageFunction.h" +#include "itkNearestNeighborInterpolateImageFunction.h" + +#include "itkImageMaskSpatialObject.h" +#include "itkRegionOfInterestImageFilter.h" +#include "itkListSample.h" + +namespace mitk +{ + + PartialVolumeAnalysisHistogramCalculator::PartialVolumeAnalysisHistogramCalculator() + : m_MaskingMode( MASKING_MODE_NONE ), + m_MaskingModeChanged( false ), + m_NumberOfBins(256), + m_UpsamplingFactor(1), + m_GaussianSigma(0) + { + m_EmptyHistogram = HistogramType::New(); + HistogramType::SizeType histogramSize; + histogramSize.Fill( 256 ); + m_EmptyHistogram->Initialize( histogramSize ); + + m_EmptyStatistics.Reset(); + } + + + PartialVolumeAnalysisHistogramCalculator::~PartialVolumeAnalysisHistogramCalculator() + { + } + + + void PartialVolumeAnalysisHistogramCalculator::SetImage( const mitk::Image *image ) + { + if ( m_Image != image ) + { + m_Image = image; + this->Modified(); + + m_ImageStatisticsTimeStamp.Modified(); + m_ImageStatisticsCalculationTriggerBool = true; + } + } + + void PartialVolumeAnalysisHistogramCalculator::AddAdditionalResamplingImage( const mitk::Image *image ) + { + m_AdditionalResamplingImages.push_back(image); + this->Modified(); + m_ImageStatisticsTimeStamp.Modified(); + m_ImageStatisticsCalculationTriggerBool = true; + } + + void PartialVolumeAnalysisHistogramCalculator::SetModified( ) + { + this->Modified(); + m_ImageStatisticsTimeStamp.Modified(); + m_ImageStatisticsCalculationTriggerBool = true; + m_MaskedImageStatisticsTimeStamp.Modified(); + m_MaskedImageStatisticsCalculationTriggerBool = true; + m_PlanarFigureStatisticsTimeStamp.Modified(); + m_PlanarFigureStatisticsCalculationTriggerBool = true; + } + + void PartialVolumeAnalysisHistogramCalculator::SetImageMask( const mitk::Image *imageMask ) + { + if ( m_Image.IsNull() ) + { + itkExceptionMacro( << "Image needs to be set first!" ); + } + + if ( m_Image->GetTimeSteps() != imageMask->GetTimeSteps() ) + { + itkExceptionMacro( << "Image and image mask need to have equal number of time steps!" ); + } + + if ( m_ImageMask != imageMask ) + { + m_ImageMask = imageMask; + this->Modified(); + + m_MaskedImageStatisticsTimeStamp.Modified(); + m_MaskedImageStatisticsCalculationTriggerBool = true; + } + } + + + void PartialVolumeAnalysisHistogramCalculator::SetPlanarFigure( const mitk::PlanarFigure *planarFigure ) + { + if ( m_Image.IsNull() ) + { + itkExceptionMacro( << "Image needs to be set first!" ); + } + + if ( m_PlanarFigure != planarFigure ) + { + m_PlanarFigure = planarFigure; + this->Modified(); + + m_PlanarFigureStatisticsTimeStamp.Modified(); + m_PlanarFigureStatisticsCalculationTriggerBool = true; + } + } + + + void PartialVolumeAnalysisHistogramCalculator::SetMaskingMode( unsigned int mode ) + { + if ( m_MaskingMode != mode ) + { + m_MaskingMode = mode; + m_MaskingModeChanged = true; + this->Modified(); + } + } + + + void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToNone() + { + if ( m_MaskingMode != MASKING_MODE_NONE ) + { + m_MaskingMode = MASKING_MODE_NONE; + m_MaskingModeChanged = true; + this->Modified(); + } + } + + + void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToImage() + { + if ( m_MaskingMode != MASKING_MODE_IMAGE ) + { + m_MaskingMode = MASKING_MODE_IMAGE; + m_MaskingModeChanged = true; + this->Modified(); + } + } + + + void PartialVolumeAnalysisHistogramCalculator::SetMaskingModeToPlanarFigure() + { + if ( m_MaskingMode != MASKING_MODE_PLANARFIGURE ) + { + m_MaskingMode = MASKING_MODE_PLANARFIGURE; + m_MaskingModeChanged = true; + this->Modified(); + } + } + + + + bool PartialVolumeAnalysisHistogramCalculator::ComputeStatistics() + { + + MITK_INFO << "ComputeStatistics() start"; + + if ( m_Image.IsNull() ) + { + itkExceptionMacro( << "Image not set!" ); + } + + if ( m_Image->GetReferenceCount() == 1 ) + { + MITK_INFO << "No Stats calculated; no one else holds a reference on it"; + return false; + } + + // If a mask was set but we are the only ones to still hold a reference on + // it, delete it. + if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() == 1) ) + { + m_ImageMask = NULL; + } + + + // Check if statistics is already up-to-date + unsigned long imageMTime = m_ImageStatisticsTimeStamp.GetMTime(); + unsigned long maskedImageMTime = m_MaskedImageStatisticsTimeStamp.GetMTime(); + unsigned long planarFigureMTime = m_PlanarFigureStatisticsTimeStamp.GetMTime(); + + bool imageStatisticsCalculationTrigger = m_ImageStatisticsCalculationTriggerBool; + bool maskedImageStatisticsCalculationTrigger = m_MaskedImageStatisticsCalculationTriggerBool; + bool planarFigureStatisticsCalculationTrigger = m_PlanarFigureStatisticsCalculationTriggerBool; + + if ( /*prevent calculation without mask*/ m_MaskingMode == MASKING_MODE_NONE || ( + ((m_MaskingMode != MASKING_MODE_NONE) || (imageMTime > m_Image->GetMTime() && !imageStatisticsCalculationTrigger)) + && ((m_MaskingMode != MASKING_MODE_IMAGE) || (m_ImageMask.IsNotNull() && maskedImageMTime > m_ImageMask->GetMTime() && !maskedImageStatisticsCalculationTrigger)) + && ((m_MaskingMode != MASKING_MODE_PLANARFIGURE) || (m_PlanarFigure.IsNotNull() && planarFigureMTime > m_PlanarFigure->GetMTime() && !planarFigureStatisticsCalculationTrigger)) ) ) + { + MITK_INFO << "Returning, statistics already up to date!"; + if ( m_MaskingModeChanged ) + { + m_MaskingModeChanged = false; + return true; + } + else + { + return false; + } + } + + // Reset state changed flag + m_MaskingModeChanged = false; + + + // Depending on masking mode, extract and/or generate the required image + // and mask data from the user input + this->ExtractImageAndMask( ); + + + Statistics *statistics; + HistogramType::ConstPointer *histogram; + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: + statistics = &m_ImageStatistics; + histogram = &m_ImageHistogram; + + m_ImageStatisticsTimeStamp.Modified(); + m_ImageStatisticsCalculationTriggerBool = false; + break; + + case MASKING_MODE_IMAGE: + statistics = &m_MaskedImageStatistics; + histogram = &m_MaskedImageHistogram; + + m_MaskedImageStatisticsTimeStamp.Modified(); + m_MaskedImageStatisticsCalculationTriggerBool = false; + break; + + case MASKING_MODE_PLANARFIGURE: + statistics = &m_PlanarFigureStatistics; + histogram = &m_PlanarFigureHistogram; + + m_PlanarFigureStatisticsTimeStamp.Modified(); + m_PlanarFigureStatisticsCalculationTriggerBool = false; + break; + } + + // Calculate statistics and histogram(s) + if ( m_InternalImage->GetDimension() == 3 ) + { + if ( m_MaskingMode == MASKING_MODE_NONE ) + { + AccessFixedDimensionByItk_2( + m_InternalImage, + InternalCalculateStatisticsUnmasked, + 3, + *statistics, + histogram ); + } + else + { + AccessFixedDimensionByItk_3( + m_InternalImage, + InternalCalculateStatisticsMasked, + 3, + m_InternalImageMask3D.GetPointer(), + *statistics, + histogram ); + } + } + else if ( m_InternalImage->GetDimension() == 2 ) + { + if ( m_MaskingMode == MASKING_MODE_NONE ) + { + AccessFixedDimensionByItk_2( + m_InternalImage, + InternalCalculateStatisticsUnmasked, + 2, + *statistics, + histogram ); + } + else + { + AccessFixedDimensionByItk_3( + m_InternalImage, + InternalCalculateStatisticsMasked, + 2, + m_InternalImageMask2D.GetPointer(), + *statistics, + histogram ); + } + } + else + { + MITK_ERROR << "ImageStatistics: Image dimension not supported!"; + } + + // Release unused image smart pointers to free memory + // m_InternalImage = mitk::Image::Pointer(); + m_InternalImageMask3D = MaskImage3DType::Pointer(); + m_InternalImageMask2D = MaskImage2DType::Pointer(); + return true; + } + + + const PartialVolumeAnalysisHistogramCalculator::HistogramType * + PartialVolumeAnalysisHistogramCalculator::GetHistogram( ) const + { + if ( m_Image.IsNull() ) + { + return NULL; + } + + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: + return m_ImageHistogram; + + case MASKING_MODE_IMAGE: + return m_MaskedImageHistogram; + + case MASKING_MODE_PLANARFIGURE: + return m_PlanarFigureHistogram; + } + } + + + const PartialVolumeAnalysisHistogramCalculator::Statistics & + PartialVolumeAnalysisHistogramCalculator::GetStatistics( ) const + { + if ( m_Image.IsNull() ) + { + return m_EmptyStatistics; + } + + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + default: + return m_ImageStatistics; + + case MASKING_MODE_IMAGE: + return m_MaskedImageStatistics; + + case MASKING_MODE_PLANARFIGURE: + return m_PlanarFigureStatistics; + } + } + + + void PartialVolumeAnalysisHistogramCalculator::ExtractImageAndMask( ) + { + + MITK_INFO << "ExtractImageAndMask( ) start"; + + if ( m_Image.IsNull() ) + { + throw std::runtime_error( "Error: image empty!" ); + } + + mitk::Image *timeSliceImage = const_cast(m_Image.GetPointer());//imageTimeSelector->GetOutput(); + + switch ( m_MaskingMode ) + { + case MASKING_MODE_NONE: + { + m_InternalImage = timeSliceImage; + int s = m_AdditionalResamplingImages.size(); + m_InternalAdditionalResamplingImages.resize(s); + for(int i=0; i(m_AdditionalResamplingImages[i].GetPointer()); + } + m_InternalImageMask2D = NULL; + m_InternalImageMask3D = NULL; + break; + } + + case MASKING_MODE_IMAGE: + { + if ( m_ImageMask.IsNotNull() && (m_ImageMask->GetReferenceCount() > 1) ) + { + ImageTimeSelector::Pointer maskedImageTimeSelector = ImageTimeSelector::New(); + maskedImageTimeSelector->SetInput( m_ImageMask ); + maskedImageTimeSelector->SetTimeNr( 0 ); + maskedImageTimeSelector->UpdateLargestPossibleRegion(); + mitk::Image *timeSliceMaskedImage = maskedImageTimeSelector->GetOutput(); + + InternalMaskImage(timeSliceMaskedImage); + + if(m_UpsamplingFactor != 1) + { + InternalResampleImage(m_InternalImageMask3D); + } + + AccessFixedDimensionByItk_1( + timeSliceImage, + InternalResampleImageFromMask, 3, -1 ); + + int s = m_AdditionalResamplingImages.size(); + m_InternalAdditionalResamplingImages.resize(s); + for(int i=0; iIsClosed() ) + { + throw std::runtime_error( "Masking not possible for non-closed figures" ); + } + + const Geometry3D *imageGeometry = timeSliceImage->GetUpdatedGeometry(); + if ( imageGeometry == NULL ) + { + throw std::runtime_error( "Image geometry invalid!" ); + } + + const Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); + if ( planarFigureGeometry2D == NULL ) + { + throw std::runtime_error( "Planar-Figure not yet initialized!" ); + } + + const PlaneGeometry *planarFigureGeometry = + dynamic_cast< const PlaneGeometry * >( planarFigureGeometry2D ); + if ( planarFigureGeometry == NULL ) + { + throw std::runtime_error( "Non-planar planar figures not supported!" ); + } + +// unsigned int axis = 2; +// unsigned int slice = 0; + + AccessFixedDimensionByItk_3( + timeSliceImage, + InternalReorientImagePlane, 3, + timeSliceImage->GetGeometry(), + m_PlanarFigure->GetGeometry(), -1 ); + + AccessFixedDimensionByItk_1( + m_InternalImage, + InternalCalculateMaskFromPlanarFigure, + 3, 2 ); + + int s = m_AdditionalResamplingImages.size(); + for(int i=0; iGetGeometry(), + m_PlanarFigure->GetGeometry(), i ); + + AccessFixedDimensionByItk_1( + m_InternalAdditionalResamplingImages[i], + InternalCropAdditionalImage, 3, i ); + } + } + } + + } + + + bool PartialVolumeAnalysisHistogramCalculator::GetPrincipalAxis( + const Geometry3D *geometry, Vector3D vector, + unsigned int &axis ) + { + vector.Normalize(); + for ( unsigned int i = 0; i < 3; ++i ) + { + Vector3D axisVector = geometry->GetAxisVector( i ); + axisVector.Normalize(); + + if ( fabs( fabs( axisVector * vector ) - 1.0) < mitk::eps ) + { + axis = i; + return true; + } + } + + return false; + } + + void PartialVolumeAnalysisHistogramCalculator::InternalMaskImage( + mitk::Image *image ) + { + + typedef itk::ImageMaskSpatialObject<3> ImageMaskSpatialObject; + typedef itk::Image< unsigned char, 3 > ImageType; + typedef itk::ImageRegion<3> RegionType; + + typedef mitk::ImageToItk CastType; + CastType::Pointer caster = CastType::New(); + caster->SetInput(image); + caster->Update(); + + ImageMaskSpatialObject::Pointer maskSO = ImageMaskSpatialObject::New(); + maskSO->SetImage ( caster->GetOutput() ); + m_InternalMask3D = + maskSO->GetAxisAlignedBoundingBoxRegion(); + + MITK_INFO << "Bounding Box Region: " << m_InternalMask3D; + + typedef itk::RegionOfInterestImageFilter< ImageType, MaskImage3DType > ROIFilterType; + ROIFilterType::Pointer roi = ROIFilterType::New(); + roi->SetRegionOfInterest(m_InternalMask3D); + roi->SetInput(caster->GetOutput()); + roi->Update(); + + m_InternalImageMask3D = roi->GetOutput(); + + MITK_INFO << "Created m_InternalImageMask3D"; + + } + + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisHistogramCalculator::InternalReorientImagePlane( + const itk::Image< TPixel, VImageDimension > *image, + mitk::Geometry3D* imggeo, mitk::Geometry3D* planegeo3D, int additionalIndex ) + { + + MITK_INFO << "InternalReorientImagePlane() start"; + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< float, VImageDimension > FloatImageType; + + typedef itk::ResampleImageFilter ResamplerType; + typename ResamplerType::Pointer resampler = ResamplerType::New(); + + mitk::PlaneGeometry* planegeo = dynamic_cast(planegeo3D); + + float upsamp = m_UpsamplingFactor; + float gausssigma = m_GaussianSigma; + + // Spacing + typename ResamplerType::SpacingType spacing = planegeo->GetSpacing(); + spacing[0] = image->GetSpacing()[0] / upsamp; + spacing[1] = image->GetSpacing()[1] / upsamp; + spacing[2] = image->GetSpacing()[2]; + resampler->SetOutputSpacing( spacing ); + + // Size + typename ResamplerType::SizeType size; + size[0] = planegeo->GetParametricExtentInMM(0) / spacing[0]; + size[1] = planegeo->GetParametricExtentInMM(1) / spacing[1]; + size[2] = 1; + resampler->SetSize( size ); + + // Origin + typename mitk::Point3D orig = planegeo->GetOrigin(); + typename mitk::Point3D corrorig; + planegeo3D->WorldToIndex(orig,corrorig); + corrorig[0] += 0.5/upsamp; + corrorig[1] += 0.5/upsamp; + corrorig[2] += 0; + planegeo3D->IndexToWorld(corrorig,corrorig); + resampler->SetOutputOrigin(corrorig ); + + // Direction + typename ResamplerType::DirectionType direction; + typename mitk::AffineTransform3D::MatrixType matrix = planegeo->GetIndexToWorldTransform()->GetMatrix(); + for(int c=0; cSetOutputDirection( direction ); + + // Gaussian interpolation + if(gausssigma != 0) + { + double sigma[3]; + for( unsigned int d = 0; d < 3; d++ ) + { + sigma[d] = gausssigma * image->GetSpacing()[d]; + } + double alpha = 2.0; + + typedef itk::GaussianInterpolateImageFunction + GaussianInterpolatorType; + + typename GaussianInterpolatorType::Pointer interpolator + = GaussianInterpolatorType::New(); + + interpolator->SetInputImage( image ); + interpolator->SetParameters( sigma, alpha ); + + resampler->SetInterpolator( interpolator ); + } + else + { + // typedef typename itk::BSplineInterpolateImageFunction + // InterpolatorType; + typedef typename itk::LinearInterpolateImageFunction InterpolatorType; + + typename InterpolatorType::Pointer interpolator + = InterpolatorType::New(); + + interpolator->SetInputImage( image ); + + resampler->SetInterpolator( interpolator ); + + } + + // Other resampling options + resampler->SetInput( image ); + resampler->SetDefaultPixelValue(0); + + MITK_INFO << "Resampling requested image plane ... "; + resampler->Update(); + MITK_INFO << " ... done"; + + if(additionalIndex < 0) + { + this->m_InternalImage = mitk::Image::New(); + this->m_InternalImage->InitializeByItk( resampler->GetOutput() ); + this->m_InternalImage->SetVolume( resampler->GetOutput()->GetBufferPointer() ); + } + else + { + unsigned int myIndex = additionalIndex; + this->m_InternalAdditionalResamplingImages.push_back(mitk::Image::New()); + this->m_InternalAdditionalResamplingImages[myIndex]->InitializeByItk( resampler->GetOutput() ); + this->m_InternalAdditionalResamplingImages[myIndex]->SetVolume( resampler->GetOutput()->GetBufferPointer() ); + } + + } + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisHistogramCalculator::InternalResampleImageFromMask( + const itk::Image< TPixel, VImageDimension > *image, int additionalIndex ) + { + typedef itk::Image< TPixel, VImageDimension > ImageType; + + typename ImageType::Pointer outImage = ImageType::New(); + outImage->SetSpacing( m_InternalImageMask3D->GetSpacing() ); // Set the image spacing + outImage->SetOrigin( m_InternalImageMask3D->GetOrigin() ); // Set the image origin + outImage->SetDirection( m_InternalImageMask3D->GetDirection() ); // Set the image direction + outImage->SetRegions( m_InternalImageMask3D->GetLargestPossibleRegion() ); + outImage->Allocate(); + outImage->FillBuffer(0); + + typedef itk::InterpolateImageFunction + BaseInterpType; + typedef itk::GaussianInterpolateImageFunction + GaussianInterpolatorType; + typedef itk::BSplineInterpolateImageFunction + BSplineInterpolatorType; + + typename BaseInterpType::Pointer interpolator; + + if(m_GaussianSigma != 0) + { + double sigma[3]; + for( unsigned int d = 0; d < 3; d++ ) + { + sigma[d] = m_GaussianSigma * image->GetSpacing()[d]; + } + double alpha = 2.0; + + interpolator = GaussianInterpolatorType::New(); + dynamic_cast(interpolator.GetPointer())->SetParameters( sigma, alpha ); + + } + else + { + interpolator = BSplineInterpolatorType::New(); + } + + interpolator->SetInputImage( image ); + + itk::ImageRegionConstIterator + itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); + itk::ImageRegionIterator + itimage(outImage, outImage->GetLargestPossibleRegion()); + + itmask = itmask.Begin(); + itimage = itimage.Begin(); + + itk::Point< double, 3 > point; + itk::ContinuousIndex< double, 3 > index; + while( !itmask.IsAtEnd() ) + { + if(itmask.Get() != 0) + { + outImage->TransformIndexToPhysicalPoint (itimage.GetIndex(), point); + image->TransformPhysicalPointToContinuousIndex(point, index); + itimage.Set(interpolator->EvaluateAtContinuousIndex(index)); + } + + ++itmask; + ++itimage; + } + + if(additionalIndex < 0) + { + this->m_InternalImage = mitk::Image::New(); + this->m_InternalImage->InitializeByItk( outImage.GetPointer() ); + this->m_InternalImage->SetVolume( outImage->GetBufferPointer() ); + } + else + { + this->m_InternalAdditionalResamplingImages[additionalIndex] = mitk::Image::New(); + this->m_InternalAdditionalResamplingImages[additionalIndex]->InitializeByItk( outImage.GetPointer() ); + this->m_InternalAdditionalResamplingImages[additionalIndex]->SetVolume( outImage->GetBufferPointer() ); + } + + } + + void PartialVolumeAnalysisHistogramCalculator::InternalResampleImage( + const MaskImage3DType *image ) + { + typedef itk::ResampleImageFilter ResamplerType; + ResamplerType::Pointer resampler = ResamplerType::New(); + + // Size + ResamplerType::SizeType size; + size[0] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[0]; + size[1] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[1]; + size[2] = m_UpsamplingFactor * image->GetLargestPossibleRegion().GetSize()[2];; + resampler->SetSize( size ); + + // Origin + mitk::Point3D orig = image->GetOrigin(); + resampler->SetOutputOrigin(orig ); + + // Spacing + ResamplerType::SpacingType spacing; + spacing[0] = image->GetSpacing()[0] / m_UpsamplingFactor; + spacing[1] = image->GetSpacing()[1] / m_UpsamplingFactor; + spacing[2] = image->GetSpacing()[2] / m_UpsamplingFactor; + resampler->SetOutputSpacing( spacing ); + + resampler->SetOutputDirection( image->GetDirection() ); + + typedef itk::NearestNeighborInterpolateImageFunction + InterpolatorType; + + InterpolatorType::Pointer interpolator + = InterpolatorType::New(); + + interpolator->SetInputImage( image ); + + resampler->SetInterpolator( interpolator ); + + // Other resampling options + resampler->SetInput( image ); + resampler->SetDefaultPixelValue(0); + resampler->Update(); + + m_InternalImageMask3D = resampler->GetOutput(); + + } + + + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisHistogramCalculator::InternalCalculateStatisticsUnmasked( + const itk::Image< TPixel, VImageDimension > *image, + Statistics &statistics, + typename HistogramType::ConstPointer *histogram ) + { + MITK_INFO << "InternalCalculateStatisticsUnmasked()"; + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned char, VImageDimension > MaskImageType; + typedef typename ImageType::IndexType IndexType; + + // Progress listening... + typedef itk::SimpleMemberCommand< PartialVolumeAnalysisHistogramCalculator > ITKCommandType; + ITKCommandType::Pointer progressListener; + progressListener = ITKCommandType::New(); + progressListener->SetCallbackFunction( this, + &PartialVolumeAnalysisHistogramCalculator::UnmaskedStatisticsProgressUpdate ); + + // Issue 100 artificial progress events since ScalarIMageToHistogramGenerator + // does not (yet?) support progress reporting + this->InvokeEvent( itk::StartEvent() ); + for ( unsigned int i = 0; i < 100; ++i ) + { + this->UnmaskedStatisticsProgressUpdate(); + } + + // Calculate statistics (separate filter) + typedef itk::StatisticsImageFilter< ImageType > StatisticsFilterType; + typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); + statisticsFilter->SetInput( image ); + unsigned long observerTag = statisticsFilter->AddObserver( + itk::ProgressEvent(), progressListener ); + + statisticsFilter->Update(); + + statisticsFilter->RemoveObserver( observerTag ); + + + this->InvokeEvent( itk::EndEvent() ); + + statistics.N = image->GetBufferedRegion().GetNumberOfPixels(); + statistics.Min = statisticsFilter->GetMinimum(); + statistics.Max = statisticsFilter->GetMaximum(); + statistics.Mean = statisticsFilter->GetMean(); + statistics.Median = 0.0; + statistics.Sigma = statisticsFilter->GetSigma(); + statistics.RMS = sqrt( statistics.Mean * statistics.Mean + + statistics.Sigma * statistics.Sigma ); + + typename ImageType::Pointer inImage = const_cast(image); + + // Calculate histogram + typedef itk::Statistics::ScalarImageToHistogramGenerator< ImageType > + HistogramGeneratorType; + typename HistogramGeneratorType::Pointer histogramGenerator = HistogramGeneratorType::New(); + histogramGenerator->SetInput( inImage ); + histogramGenerator->SetMarginalScale( 10 ); // Defines y-margin width of histogram + histogramGenerator->SetNumberOfBins( m_NumberOfBins ); // CT range [-1024, +2048] --> bin size 4 values + histogramGenerator->SetHistogramMin( statistics.Min ); + histogramGenerator->SetHistogramMax( statistics.Max ); + histogramGenerator->Compute(); + *histogram = histogramGenerator->GetOutput(); + } + + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisHistogramCalculator::InternalCalculateStatisticsMasked( + const itk::Image< TPixel, VImageDimension > *image, + itk::Image< unsigned char, VImageDimension > *maskImage, + Statistics &statistics, + typename HistogramType::ConstPointer *histogram ) + { + MITK_INFO << "InternalCalculateStatisticsMasked() start"; + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::Image< unsigned char, VImageDimension > MaskImageType; + typedef typename ImageType::IndexType IndexType; + + // generate a list sample of angles at positions + // where the fiber-prob is higher than .2*maxprob + typedef TPixel MeasurementType; + const unsigned int MeasurementVectorLength = 1; + typedef itk::Vector< MeasurementType , MeasurementVectorLength > + MeasurementVectorType; + typedef itk::Statistics::ListSample< MeasurementVectorType > ListSampleType; + typename ListSampleType::Pointer listSample = ListSampleType::New(); + listSample->SetMeasurementVectorSize( MeasurementVectorLength ); + + itk::ImageRegionConstIterator + itmask(maskImage, maskImage->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator + itimage(image, image->GetLargestPossibleRegion()); + + itmask = itmask.Begin(); + itimage = itimage.Begin(); + + while( !itmask.IsAtEnd() ) + { + if(itmask.Get() != 0) + { + // apend to list + MeasurementVectorType mv; + mv[0] = ( MeasurementType ) itimage.Get(); + listSample->PushBack(mv); + } + ++itmask; + ++itimage; + } + + // generate a histogram from the list sample + typedef float HistogramMeasurementType; + typedef itk::Statistics::ListSampleToHistogramGenerator + < ListSampleType, HistogramMeasurementType, + itk::Statistics::DenseFrequencyContainer, + MeasurementVectorLength > GeneratorType; + typename GeneratorType::Pointer generator = GeneratorType::New(); + typename GeneratorType::HistogramType::SizeType size; + size.Fill(m_NumberOfBins); + generator->SetNumberOfBins( size ); + generator->SetListSample( listSample ); + generator->SetMarginalScale( 10.0 ); + generator->Update(); + *histogram = generator->GetOutput(); + + } + + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisHistogramCalculator::InternalCropAdditionalImage( + itk::Image< TPixel, VImageDimension > *image, int additionalIndex ) + { + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > ROIFilterType; + typename ROIFilterType::Pointer roi = ROIFilterType::New(); + roi->SetRegionOfInterest(m_CropRegion); + roi->SetInput(image); + roi->Update(); + + m_InternalAdditionalResamplingImages[additionalIndex] = mitk::Image::New(); + m_InternalAdditionalResamplingImages[additionalIndex]->InitializeByItk(roi->GetOutput()); + m_InternalAdditionalResamplingImages[additionalIndex]->SetVolume(roi->GetOutput()->GetBufferPointer()); + } + + template < typename TPixel, unsigned int VImageDimension > + void PartialVolumeAnalysisHistogramCalculator::InternalCalculateMaskFromPlanarFigure( + itk::Image< TPixel, VImageDimension > *image, unsigned int axis ) + { + + MITK_INFO << "InternalCalculateMaskFromPlanarFigure() start"; + + typedef itk::Image< TPixel, VImageDimension > ImageType; + typedef itk::CastImageFilter< ImageType, MaskImage3DType > CastFilterType; + + // Generate mask image as new image with same header as input image and + // initialize with "1". + MaskImage3DType::Pointer newMaskImage = MaskImage3DType::New(); + newMaskImage->SetSpacing( image->GetSpacing() ); // Set the image spacing + newMaskImage->SetOrigin( image->GetOrigin() ); // Set the image origin + newMaskImage->SetDirection( image->GetDirection() ); // Set the image direction + newMaskImage->SetRegions( image->GetLargestPossibleRegion() ); + newMaskImage->Allocate(); + newMaskImage->FillBuffer( 1 ); + + // Generate VTK polygon from (closed) PlanarFigure polyline + // (The polyline points are shifted by -0.5 in z-direction to make sure + // that the extrusion filter, which afterwards elevates all points by +0.5 + // in z-direction, creates a 3D object which is cut by the the plane z=0) + const mitk::Geometry2D *planarFigureGeometry2D = m_PlanarFigure->GetGeometry2D(); + const typename PlanarFigure::PolyLineType planarFigurePolyline = m_PlanarFigure->GetPolyLine( 0 ); + const mitk::Geometry3D *imageGeometry3D = m_InternalImage->GetGeometry( 0 ); + + vtkPolyData *polyline = vtkPolyData::New(); + polyline->Allocate( 1, 1 ); + + // Determine x- and y-dimensions depending on principal axis + int i0, i1; + switch ( axis ) + { + case 0: + i0 = 1; + i1 = 2; + break; + + case 1: + i0 = 0; + i1 = 2; + break; + + case 2: + default: + i0 = 0; + i1 = 1; + break; + } + + // Create VTK polydata object of polyline contour + bool outOfBounds = false; + vtkPoints *points = vtkPoints::New(); + typename PlanarFigure::PolyLineType::const_iterator it; + for ( it = planarFigurePolyline.begin(); + it != planarFigurePolyline.end(); + ++it ) + { + Point3D point3D; + + // Convert 2D point back to the local index coordinates of the selected + // image + mitk::Point2D point2D = it->Point; + planarFigureGeometry2D->WorldToIndex(point2D, point2D); + point2D[0] -= 0.5/m_UpsamplingFactor; + point2D[1] -= 0.5/m_UpsamplingFactor; + planarFigureGeometry2D->IndexToWorld(point2D, point2D); + planarFigureGeometry2D->Map( point2D, point3D ); + + // Polygons (partially) outside of the image bounds can not be processed + // further due to a bug in vtkPolyDataToImageStencil + if ( !imageGeometry3D->IsInside( point3D ) ) + { + outOfBounds = true; + } + + imageGeometry3D->WorldToIndex( point3D, point3D ); + point3D[i0] += 0.5; + point3D[i1] += 0.5; + + // Add point to polyline array + points->InsertNextPoint( point3D[i0], point3D[i1], -0.5 ); + } + polyline->SetPoints( points ); + points->Delete(); + + if ( outOfBounds ) + { + polyline->Delete(); + throw std::runtime_error( "Figure at least partially outside of image bounds!" ); + } + + unsigned int numberOfPoints = planarFigurePolyline.size(); + vtkIdType *ptIds = new vtkIdType[numberOfPoints]; + for ( vtkIdType i = 0; i < numberOfPoints; ++i ) + { + ptIds[i] = i; + } + polyline->InsertNextCell( VTK_POLY_LINE, numberOfPoints, ptIds ); + + + // Extrude the generated contour polygon + vtkLinearExtrusionFilter *extrudeFilter = vtkLinearExtrusionFilter::New(); + extrudeFilter->SetInput( polyline ); + extrudeFilter->SetScaleFactor( 1 ); + extrudeFilter->SetExtrusionTypeToNormalExtrusion(); + extrudeFilter->SetVector( 0.0, 0.0, 1.0 ); + + // Make a stencil from the extruded polygon + vtkPolyDataToImageStencil *polyDataToImageStencil = vtkPolyDataToImageStencil::New(); + polyDataToImageStencil->SetInput( extrudeFilter->GetOutput() ); + + + + // Export from ITK to VTK (to use a VTK filter) + typedef itk::VTKImageImport< MaskImage3DType > ImageImportType; + typedef itk::VTKImageExport< MaskImage3DType > ImageExportType; + + typename ImageExportType::Pointer itkExporter = ImageExportType::New(); + itkExporter->SetInput( newMaskImage ); + + vtkImageImport *vtkImporter = vtkImageImport::New(); + this->ConnectPipelines( itkExporter, vtkImporter ); + vtkImporter->Update(); + + + // Apply the generated image stencil to the input image + vtkImageStencil *imageStencilFilter = vtkImageStencil::New(); + imageStencilFilter->SetInput( vtkImporter->GetOutput() ); + imageStencilFilter->SetStencil( polyDataToImageStencil->GetOutput() ); + imageStencilFilter->ReverseStencilOff(); + imageStencilFilter->SetBackgroundValue( 0 ); + imageStencilFilter->Update(); + + + // Export from VTK back to ITK + vtkImageExport *vtkExporter = vtkImageExport::New(); + vtkExporter->SetInput( imageStencilFilter->GetOutput() ); + vtkExporter->Update(); + + typename ImageImportType::Pointer itkImporter = ImageImportType::New(); + this->ConnectPipelines( vtkExporter, itkImporter ); + itkImporter->Update(); + + // calculate cropping bounding box + m_InternalImageMask3D = itkImporter->GetOutput(); + m_InternalImageMask3D->SetDirection(image->GetDirection()); + + itk::ImageRegionConstIterator + itmask(m_InternalImageMask3D, m_InternalImageMask3D->GetLargestPossibleRegion()); + itk::ImageRegionIterator + itimage(image, image->GetLargestPossibleRegion()); + + itmask = itmask.Begin(); + itimage = itimage.Begin(); + + typename ImageType::SizeType lowersize = {{9999999999.0,9999999999.0,9999999999.0}}; + typename ImageType::SizeType uppersize = {{0,0,0}}; + while( !itmask.IsAtEnd() ) + { + if(itmask.Get() == 0) + { + itimage.Set(0); + } + else + { + typename ImageType::IndexType index = itimage.GetIndex(); + typename ImageType::SizeType signedindex; + signedindex[0] = index[0]; + signedindex[1] = index[1]; + signedindex[2] = index[2]; + + lowersize[0] = signedindex[0] < lowersize[0] ? signedindex[0] : lowersize[0]; + lowersize[1] = signedindex[1] < lowersize[1] ? signedindex[1] : lowersize[1]; + lowersize[2] = signedindex[2] < lowersize[2] ? signedindex[2] : lowersize[2]; + + uppersize[0] = signedindex[0] > uppersize[0] ? signedindex[0] : uppersize[0]; + uppersize[1] = signedindex[1] > uppersize[1] ? signedindex[1] : uppersize[1]; + uppersize[2] = signedindex[2] > uppersize[2] ? signedindex[2] : uppersize[2]; + } + + ++itmask; + ++itimage; + } + + typename ImageType::IndexType index; + index[0] = lowersize[0]; + index[1] = lowersize[1]; + index[2] = lowersize[2]; + + typename ImageType::SizeType size; + size[0] = uppersize[0] - lowersize[0] + 1; + size[1] = uppersize[1] - lowersize[1] + 1; + size[2] = uppersize[2] - lowersize[2] + 1; + + m_CropRegion = itk::ImageRegion<3>(index, size); + + // crop internal image + typedef itk::RegionOfInterestImageFilter< ImageType, ImageType > ROIFilterType; + typename ROIFilterType::Pointer roi = ROIFilterType::New(); + roi->SetRegionOfInterest(m_CropRegion); + roi->SetInput(image); + roi->Update(); + + m_InternalImage = mitk::Image::New(); + m_InternalImage->InitializeByItk(roi->GetOutput()); + m_InternalImage->SetVolume(roi->GetOutput()->GetBufferPointer()); + + // crop internal mask + typedef itk::RegionOfInterestImageFilter< MaskImage3DType, MaskImage3DType > ROIMaskFilterType; + typename ROIMaskFilterType::Pointer roi2 = ROIMaskFilterType::New(); + roi2->SetRegionOfInterest(m_CropRegion); + roi2->SetInput(m_InternalImageMask3D); + roi2->Update(); + m_InternalImageMask3D = roi2->GetOutput(); + + // Clean up VTK objects + polyline->Delete(); + extrudeFilter->Delete(); + polyDataToImageStencil->Delete(); + vtkImporter->Delete(); + imageStencilFilter->Delete(); + //vtkExporter->Delete(); // TODO: crashes when outcommented; memory leak?? + delete[] ptIds; + + } + + + void PartialVolumeAnalysisHistogramCalculator::UnmaskedStatisticsProgressUpdate() + { + // Need to throw away every second progress event to reach a final count of + // 100 since two consecutive filters are used in this case + static int updateCounter = 0; + if ( updateCounter++ % 2 == 0 ) + { + this->InvokeEvent( itk::ProgressEvent() ); + } + } + + + void PartialVolumeAnalysisHistogramCalculator::MaskedStatisticsProgressUpdate() + { + this->InvokeEvent( itk::ProgressEvent() ); + } + + +} diff --git a/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h new file mode 100644 index 0000000000..92fb3f6a41 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h @@ -0,0 +1,353 @@ +/*======================================================================== +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + + +#ifndef _MITK_PartialVolumeAnalysisHistogramCalculator_H +#define _MITK_PartialVolumeAnalysisHistogramCalculator_H + +#include "MitkDiffusionImagingExports.h" + +#include +#include +#include + + +#include "mitkImage.h" +#include "mitkImageTimeSelector.h" +#include "mitkPlanarFigure.h" + + + +namespace mitk +{ + + /** + * \brief Class for calculating statistics and histogram for an (optionally + * masked) image. + * + * Images can be masked by either a (binary) image (of the same dimensions as + * the original image) or by a closed mitk::PlanarFigure, e.g. a circle or + * polygon. When masking with a planar figure, the slice corresponding to the + * plane containing the figure is extracted and then clipped with contour + * defined by the figure. Planar figures need to be aligned along the main axes + * of the image (transversal, sagittal, coronal). Planar figures on arbitrary + * rotated planes are not supported. + * + * For each operating mode (no masking, masking by image, masking by planar + * figure), the calculated statistics and histogram are cached so that, when + * switching back and forth between operation modes without modifying mask or + * image, the information doesn't need to be recalculated. + * + * Note: currently time-resolved and multi-channel pictures are not properly + * supported. + */ + class MitkDiffusionImaging_EXPORT PartialVolumeAnalysisHistogramCalculator : public itk::Object + { + public: + + enum + { + MASKING_MODE_NONE = 0, + MASKING_MODE_IMAGE, + MASKING_MODE_PLANARFIGURE + }; + + typedef mitk::Image::HistogramType HistogramType; + typedef mitk::Image::HistogramType::ConstIterator HistogramConstIteratorType; + + struct Statistics + { + unsigned int N; + double Min; + double Max; + double Mean; + double Median; + double Variance; + double Sigma; + double RMS; + + void Reset() + { + N = 0; + Min = 0.0; + Max = 0.0; + Mean = 0.0; + Median = 0.0; + Variance = 0.0; + Sigma = 0.0; + RMS = 0.0; + } + }; + + typedef Statistics StatisticsType; + + typedef itk::TimeStamp TimeStampType; + typedef bool BoolType; + + typedef itk::Image< unsigned char, 3 > MaskImage3DType; + typedef itk::Image< unsigned char, 2 > MaskImage2DType; + + typedef itk::Image< float, 2 > InternalImage2DType; + + mitkClassMacro( PartialVolumeAnalysisHistogramCalculator, itk::Object ) + itkNewMacro( PartialVolumeAnalysisHistogramCalculator ) + + /** \brief Set image from which to compute statistics. */ + void SetImage( const mitk::Image *image ); + + /** \brief Set binary image for masking. */ + void SetImageMask( const mitk::Image *imageMask ); + + /** \brief Set planar figure for masking. */ + void SetPlanarFigure( const mitk::PlanarFigure *planarFigure ); + + /** \brief Set image for which the same resampling will be applied. + and available via GetAdditionalResampledImage() */ + void AddAdditionalResamplingImage( const mitk::Image *image ); + + /** \brief Set/Get operation mode for masking */ + void SetMaskingMode( unsigned int mode ); + + /** \brief Set/Get operation mode for masking */ + itkGetMacro( MaskingMode, unsigned int ); + + /** \brief Set/Get operation mode for masking */ + void SetMaskingModeToNone(); + + /** \brief Set/Get operation mode for masking */ + void SetMaskingModeToImage(); + + /** \brief Set/Get operation mode for masking */ + void SetMaskingModeToPlanarFigure(); + + /** \brief Set histogram number of bins. */ + void SetNumberOfBins( unsigned int number ) + { + if(m_NumberOfBins != number) + { + m_NumberOfBins = number; + SetModified(); + } + } + + /** \brief Get histogram number of bins. */ + unsigned int GetNumberOfBins( ) + { return m_NumberOfBins; } + + /** \brief Set upsampling factor. */ + void SetUpsamplingFactor( float number ) + { + if(m_UpsamplingFactor != number) + { + m_UpsamplingFactor = number; + SetModified(); + } + } + + /** \brief Get upsampling factor. */ + float GetUpsamplingFactor( ) + { return m_UpsamplingFactor; } + + /** \brief Set gaussian sigma. */ + void SetGaussianSigma( float number ) + { + if(m_GaussianSigma != number) + { + m_GaussianSigma = number; + SetModified(); + } + } + + /** \brief Get histogram number of bins. */ + float GetGaussianSigma( ) + { return m_GaussianSigma; } + + void SetModified(); + + /** \brief Compute statistics (together with histogram) for the current + * masking mode. + * + * Computation is not executed if statistics is already up to date. In this + * case, false is returned; otherwise, true.*/ + virtual bool ComputeStatistics( ); + + + /** \brief Retrieve the histogram depending on the current masking mode. */ + const HistogramType *GetHistogram( ) const; + + /** \brief Retrieve statistics depending on the current masking mode. */ + const Statistics &GetStatistics( ) const; + + const Image::Pointer GetInternalImage() + { + return m_InternalImage; + } + + const Image::Pointer GetInternalAdditionalResampledImage(unsigned int i) + { + if(i < m_InternalAdditionalResamplingImages.size()) + { + return m_InternalAdditionalResamplingImages[i]; + } + else + { + return NULL; + } + } + + protected: + + PartialVolumeAnalysisHistogramCalculator(); + + virtual ~PartialVolumeAnalysisHistogramCalculator(); + + /** \brief Depending on the masking mode, the image and mask from which to + * calculate statistics is extracted from the original input image and mask + * data. + * + * For example, a when using a PlanarFigure as mask, the 2D image slice + * corresponding to the PlanarFigure will be extracted from the original + * image. If masking is disabled, the original image is simply passed + * through. */ + void ExtractImageAndMask( ); + + + /** \brief If the passed vector matches any of the three principal axes + * of the passed geometry, the ínteger value corresponding to the axis + * is set and true is returned. */ + bool GetPrincipalAxis( const Geometry3D *geometry, Vector3D vector, + unsigned int &axis ); + + template < typename TPixel, unsigned int VImageDimension > + void InternalCalculateStatisticsUnmasked( + const itk::Image< TPixel, VImageDimension > *image, + Statistics &statistics, + typename HistogramType::ConstPointer *histogram ); + + template < typename TPixel, unsigned int VImageDimension > + void InternalCalculateStatisticsMasked( + const itk::Image< TPixel, VImageDimension > *image, + itk::Image< unsigned char, VImageDimension > *maskImage, + Statistics &statistics, + typename HistogramType::ConstPointer *histogram ); + + template < typename TPixel, unsigned int VImageDimension > + void InternalCalculateMaskFromPlanarFigure( + itk::Image< TPixel, VImageDimension > *image, unsigned int axis ); + + template < typename TPixel, unsigned int VImageDimension > + void InternalReorientImagePlane( + const itk::Image< TPixel, VImageDimension > *image, mitk::Geometry3D* imggeo, mitk::Geometry3D* planegeo3D, int additionalIndex ); + + template < typename TPixel, unsigned int VImageDimension > + void InternalResampleImageFromMask( + const itk::Image< TPixel, VImageDimension > *image, int additionalIndex ); + + void InternalResampleImage( + const MaskImage3DType *image/*, mitk::Geometry3D* imggeo*/ ); + + template < typename TPixel, unsigned int VImageDimension > + void InternalCropAdditionalImage( + itk::Image< TPixel, VImageDimension > *image, int additionalIndex ); + + void InternalMaskImage( mitk::Image *image ); + + /** Connection from ITK to VTK */ + template + void ConnectPipelines(ITK_Exporter exporter, VTK_Importer* importer) + { + importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); + + importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); + importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); + importer->SetSpacingCallback(exporter->GetSpacingCallback()); + importer->SetOriginCallback(exporter->GetOriginCallback()); + importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); + + importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); + + importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); + importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); + importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); + importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); + importer->SetCallbackUserData(exporter->GetCallbackUserData()); + } + + /** Connection from VTK to ITK */ + template + void ConnectPipelines(VTK_Exporter* exporter, ITK_Importer importer) + { + importer->SetUpdateInformationCallback(exporter->GetUpdateInformationCallback()); + + importer->SetPipelineModifiedCallback(exporter->GetPipelineModifiedCallback()); + importer->SetWholeExtentCallback(exporter->GetWholeExtentCallback()); + importer->SetSpacingCallback(exporter->GetSpacingCallback()); + importer->SetOriginCallback(exporter->GetOriginCallback()); + importer->SetScalarTypeCallback(exporter->GetScalarTypeCallback()); + + importer->SetNumberOfComponentsCallback(exporter->GetNumberOfComponentsCallback()); + + importer->SetPropagateUpdateExtentCallback(exporter->GetPropagateUpdateExtentCallback()); + importer->SetUpdateDataCallback(exporter->GetUpdateDataCallback()); + importer->SetDataExtentCallback(exporter->GetDataExtentCallback()); + importer->SetBufferPointerCallback(exporter->GetBufferPointerCallback()); + importer->SetCallbackUserData(exporter->GetCallbackUserData()); + } + + void UnmaskedStatisticsProgressUpdate(); + void MaskedStatisticsProgressUpdate(); + + mitk::Image::ConstPointer m_Image; + mitk::Image::ConstPointer m_ImageMask; + mitk::PlanarFigure::ConstPointer m_PlanarFigure; + + HistogramType::ConstPointer m_ImageHistogram; + HistogramType::ConstPointer m_MaskedImageHistogram; + HistogramType::ConstPointer m_PlanarFigureHistogram; + + HistogramType::Pointer m_EmptyHistogram; + + StatisticsType m_ImageStatistics; + StatisticsType m_MaskedImageStatistics; + StatisticsType m_PlanarFigureStatistics; + + Statistics m_EmptyStatistics; + + unsigned int m_MaskingMode; + bool m_MaskingModeChanged; + + mitk::Image::Pointer m_InternalImage; + MaskImage3DType::Pointer m_InternalImageMask3D; + MaskImage2DType::Pointer m_InternalImageMask2D; + itk::ImageRegion<3> m_InternalMask3D; + std::vector m_AdditionalResamplingImages; + std::vector m_InternalAdditionalResamplingImages; + + TimeStampType m_ImageStatisticsTimeStamp; + TimeStampType m_MaskedImageStatisticsTimeStamp; + TimeStampType m_PlanarFigureStatisticsTimeStamp; + + BoolType m_ImageStatisticsCalculationTriggerBool; + BoolType m_MaskedImageStatisticsCalculationTriggerBool; + BoolType m_PlanarFigureStatisticsCalculationTriggerBool; + + unsigned int m_NumberOfBins; + + float m_UpsamplingFactor; + + float m_GaussianSigma; + + itk::ImageRegion<3> m_CropRegion; + }; + +} + +#endif // #define _MITK_PartialVolumeAnalysisHistogramCalculator_H diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 817c2ad931..1fb7dd7b1c 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,143 +1,150 @@ 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 # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp - # DataStructures -> Tbss - IODataStructures/TbssImages/mitkTbssImageSource.cpp - IODataStructures/TbssImages/mitkNrrdTbssImageReader.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/mitkFiberBundleMapper3D.cpp Rendering/mitkFiberBundleXMapper3D.cpp - + # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms + Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp + Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp ) SET(H_FILES # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleMapper3D.h Rendering/mitkFiberBundleXMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/FiberBundle/itkSlowPolyLineParametricPath.h 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 - + # Tractography Tractography/itkGlobalTractographyFilter.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkTractsToProbabilityImageFilter.h 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 )