diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake index 525cd9521b..91f43c4ac5 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/files.cmake @@ -1,102 +1,111 @@ SET(SRC_CPP_FILES QmitkODFDetailsWidget.cpp QmitkODFRenderWidget.cpp QmitkPartialVolumeAnalysisWidget.cpp QmitkIVIMWidget.cpp + QmitkTbssRoiAnalysisWidget.cpp ) SET(INTERNAL_CPP_FILES mitkPluginActivator.cpp QmitkQBallReconstructionView.cpp QmitkPreprocessingView.cpp QmitkDiffusionDicomImportView.cpp QmitkDiffusionQuantificationView.cpp QmitkTensorReconstructionView.cpp QmitkDiffusionImagingPublicPerspective.cpp QmitkControlVisualizationPropertiesView.cpp QmitkODFDetailsView.cpp QmitkGibbsTrackingView.cpp QmitkStochasticFiberTrackingView.cpp QmitkFiberProcessingView.cpp QmitkFiberBundleDeveloperView.cpp QmitkPartialVolumeAnalysisView.cpp QmitkIVIMView.cpp QmitkScreenshotMaker.cpp + QmitkTractbasedSpatialStatisticsView.cpp + QmitkTbssTableModel.cpp + QmitkTbssMetaTableModel.cpp ) SET(UI_FILES src/internal/QmitkQBallReconstructionViewControls.ui src/internal/QmitkPreprocessingViewControls.ui src/internal/QmitkDiffusionDicomImportViewControls.ui src/internal/QmitkDiffusionQuantificationViewControls.ui src/internal/QmitkTensorReconstructionViewControls.ui src/internal/QmitkControlVisualizationPropertiesViewControls.ui src/internal/QmitkODFDetailsViewControls.ui src/internal/QmitkGibbsTrackingViewControls.ui src/internal/QmitkStochasticFiberTrackingViewControls.ui src/internal/QmitkFiberProcessingViewControls.ui src/internal/QmitkFiberBundleDeveloperViewControls.ui src/internal/QmitkPartialVolumeAnalysisViewControls.ui src/internal/QmitkIVIMViewControls.ui src/internal/QmitkScreenshotMakerControls.ui + src/internal/QmitkTractbasedSpatialStatisticsViewControls.ui ) SET(MOC_H_FILES src/internal/mitkPluginActivator.h src/internal/QmitkQBallReconstructionView.h src/internal/QmitkPreprocessingView.h src/internal/QmitkDiffusionDicomImportView.h src/internal/QmitkDiffusionImagingPublicPerspective.h src/internal/QmitkDiffusionQuantificationView.h src/internal/QmitkTensorReconstructionView.h src/internal/QmitkControlVisualizationPropertiesView.h src/internal/QmitkODFDetailsView.h src/QmitkODFRenderWidget.h src/QmitkODFDetailsWidget.h src/internal/QmitkGibbsTrackingView.h src/internal/QmitkStochasticFiberTrackingView.h src/internal/QmitkFiberProcessingView.h src/internal/QmitkFiberBundleDeveloperView.h src/internal/QmitkPartialVolumeAnalysisView.h src/QmitkPartialVolumeAnalysisWidget.h src/internal/QmitkIVIMView.h src/internal/QmitkScreenshotMaker.h + src/internal/QmitkTractbasedSpatialStatisticsView.h + src/QmitkTbssRoiAnalysisWidget.h ) SET(CACHED_RESOURCE_FILES # list of resource files which can be used by the plug-in # system without loading the plug-ins shared library, # for example the icon used in the menu and tabs for the # plug-in views in the workbench plugin.xml resources/preprocessing.png resources/dwiimport.png resources/quantification.png resources/reconodf.png resources/recontensor.png resources/vizControls.png resources/OdfDetails.png resources/GibbsTracking.png resources/FiberBundleOperations.png resources/PartialVolumeAnalysis_24.png resources/IVIM_48.png resources/screenshot_maker.png resources/stochFB.png + resources/tbss.png ) SET(QRC_FILES # uncomment the following line if you want to use Qt resources resources/QmitkDiffusionImaging.qrc + #resources/QmitkTractbasedSpatialStatisticsView.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 cc0063c353..466160f756 100644 --- a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/plugin.xml @@ -1,110 +1,117 @@ + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/tbss.png b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/tbss.png new file mode 100755 index 0000000000..15950524ae Binary files /dev/null and b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/resources/tbss.png differ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp new file mode 100644 index 0000000000..ba87ece4fd --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp @@ -0,0 +1,256 @@ +/*========================================================================= +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 "QmitkTbssRoiAnalysisWidget.h" + + +#include +#include +#include + +#include +#include +#include +#include +#include + + +#include +#include + + +QmitkTbssRoiAnalysisWidget::QmitkTbssRoiAnalysisWidget( QWidget * parent ) + : QmitkPlotWidget(parent) +{ + m_PlotPicker = new QwtPlotPicker(m_Plot->canvas()); + m_PlotPicker->setSelectionFlags(QwtPicker::PointSelection | QwtPicker::ClickSelection | QwtPicker::DragSelection); + m_PlotPicker->setTrackerMode(QwtPicker::ActiveOnly); +} + +std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfiles(std::string preprocessed) +{ + MITK_INFO << "make profiles!"; + std::vector< std::vector > profiles; + + //No results were preprocessed, so they must be calculated now. + if(preprocessed == "") + { + // Iterate through the 4th dim (corresponding to subjects) + // and create a profile for every subject + int size = m_Projections->GetVectorLength(); + for(int s=0; s profile; + RoiType::iterator it; + it = m_Roi.begin(); + while(it != m_Roi.end()) + { + itk::Index<3> ix = *it; + profile.push_back(m_Projections->GetPixel(ix).GetElement(s)); + it++; + } + int pSize = profile.size(); + profiles.push_back(profile); + } + } + else{ // Use preprocessed results + std::ifstream file(preprocessed.c_str()); + if(file.is_open()) + { + std::string line; + while(getline(file,line)) + { + std::vector tokens; + Tokenize(line, tokens); + std::vector::iterator it; + it = tokens.begin(); + std::vector< double > profile; + while(it != tokens.end()) + { + std::string s = *it; + profile.push_back (atof( s.c_str() ) ); + ++it; + } + profiles.push_back(profile); + } + } + } + // Calculate the averages + + // Here a check could be build in to check whether all profiles have + // the same length, but this should normally be the case if the input + // data were corrected with the TBSS Module. + + std::vector< std::vector > groupProfiles; + + std::vector< std::pair >::iterator it; + it = m_Groups.begin(); + + int c = 0; //the current profile number + + int nprof = profiles.size(); + + while(it != m_Groups.end() && profiles.size() > 0) + { + std::pair p = *it; + int size = p.second; + + //initialize a vector of the right length with zeroes + std::vector averageProfile; + for(int i=0; iClear(); + + + + m_Vals.clear(); + + std::vector v1; + + + std::vector > groupProfiles = CalculateGroupProfiles(preprocessed); + + + std::vector xAxis; + for(int i=0; iSetPlotTitle( title.c_str() ); + QPen pen( Qt::SolidLine ); + pen.setWidth(2); + + + + std::vector< std::pair >::iterator it; + it = m_Groups.begin(); + + int c = 0; //the current profile number + + QColor colors[4] = {Qt::green, Qt::blue, Qt::yellow, Qt::red}; + + while(it != m_Groups.end() && groupProfiles.size() > 0) + { + + std::pair< std::string, int > group = *it; + + pen.setColor(colors[c]); + int curveId = this->InsertCurve( group.first.c_str() ); + this->SetCurveData( curveId, xAxis, groupProfiles.at(c) ); + + + this->SetCurvePen( curveId, pen ); + + c++; + it++; + + } + + + QwtLegend *legend = new QwtLegend; + this->SetLegend(legend, QwtPlot::RightLegend, 0.5); + + + std::cout << m_Measure << std::endl; + this->m_Plot->setAxisTitle(0, m_Measure.c_str()); + this->m_Plot->setAxisTitle(3, "Position"); + + this->Replot(); + +} + +void QmitkTbssRoiAnalysisWidget::Boxplots() +{ + this->Clear(); +} + +void QmitkTbssRoiAnalysisWidget::drawBar(int x) +{ + + m_Plot->detachItems(QwtPlotItem::Rtti_PlotMarker, true); + + + QwtPlotMarker *mX = new QwtPlotMarker(); + //mX->setLabel(QString::fromLatin1("selected point")); + mX->setLabelAlignment(Qt::AlignLeft | Qt::AlignBottom); + mX->setLabelOrientation(Qt::Vertical); + mX->setLineStyle(QwtPlotMarker::VLine); + mX->setLinePen(QPen(Qt::black, 0, Qt::SolidLine)); + mX->setXValue(x); + mX->attach(m_Plot); + + + this->Replot(); + +} + +QmitkTbssRoiAnalysisWidget::~QmitkTbssRoiAnalysisWidget() +{ + delete m_PlotPicker; + +} + + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h new file mode 100644 index 0000000000..c68f4aaa50 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h @@ -0,0 +1,148 @@ +/*========================================================================= + +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 QmitkTbssRoiAnalysisWidget_H_ +#define QmitkTbssRoiAnalysisWidget_H_ + +#include "QmitkPlotWidget.h" + +#include + +//#include "QmitkHistogram.h" +#include "QmitkExtExports.h" +#include "mitkImage.h" +#include "mitkPlanarFigure.h" +#include "itkVectorImage.h" + + +//#include + +//#include +#include +#include +#include +#include +#include + +typedef itk::VectorImage VectorImageType; + +typedef std::vector< itk::Index<3> > RoiType; + +/** + * \brief Widget for displaying boxplots + * framework + */ +class DIFFUSIONIMAGING_EXPORT QmitkTbssRoiAnalysisWidget : public QmitkPlotWidget +{ + +Q_OBJECT + +public: + + + QmitkTbssRoiAnalysisWidget( QWidget * parent); + virtual ~QmitkTbssRoiAnalysisWidget(); + + void SetGroups(std::vector< std::pair > groups) + { + m_Groups = groups; + } + + void DrawProfiles(std::string preprocessed); + + void Boxplots(); + + void SetProjections(VectorImageType::Pointer projections) + { + m_Projections = projections; + } + + void SetRoi(RoiType roi) + { + m_Roi = roi; + } + + void SetStructure(std::string structure) + { + m_Structure = structure; + } + + void SetMeasure(std::string measure) + { + m_Measure = measure; + } + + QwtPlot* GetPlot() + { + return m_Plot; + } + + QwtPlotPicker* m_PlotPicker; + + + + void drawBar(int x); + + std::vector > GetVals() + { + return m_Vals; + } + + +protected: + + std::vector< std::vector > m_Vals; + + + + + std::vector< std::vector > CalculateGroupProfiles(std::string preprocessed); + + + + void Tokenize(const std::string& str, + std::vector& tokens, + const std::string& delimiters = " ") + { + // Skip delimiters at beginning. + std::string::size_type lastPos = str.find_first_not_of(delimiters, 0); + // Find first "non-delimiter". + std::string::size_type pos = str.find_first_of(delimiters, lastPos); + + while (std::string::npos != pos || std::string::npos != lastPos) + { + // Found a token, add it to the vector. + tokens.push_back(str.substr(lastPos, pos - lastPos)); + // Skip delimiters. Note the "not_of" + lastPos = str.find_first_not_of(delimiters, pos); + // Find next "non-delimiter" + pos = str.find_first_of(delimiters, lastPos); + } + } + + std::vector< std::pair > m_Groups; + + VectorImageType::Pointer m_Projections; + RoiType m_Roi; + std::string m_Structure; + std::string m_Measure; + + + +}; + +#endif diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.cpp new file mode 100644 index 0000000000..dafdc3f1f9 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.cpp @@ -0,0 +1,150 @@ + + +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +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. + +=========================================================================*/ + + + +#include "QmitkTbssMetaTableModel.h" + +QmitkTbssMetaTableModel::QmitkTbssMetaTableModel(QObject *parent) + : QAbstractTableModel(parent) +{ +} + +QmitkTbssMetaTableModel::QmitkTbssMetaTableModel(QList< QPair > pairs, QObject *parent) + : QAbstractTableModel(parent) +{ + listOfPairs=pairs; +} + +int QmitkTbssMetaTableModel::rowCount(const QModelIndex &parent) const +{ + Q_UNUSED(parent); + return listOfPairs.size(); +} + +int QmitkTbssMetaTableModel::columnCount(const QModelIndex &parent) const +{ + Q_UNUSED(parent); + return 2; +} + +QVariant QmitkTbssMetaTableModel::data(const QModelIndex &index, int role) const +{ + if (!index.isValid()) + return QVariant(); + + if (index.row() >= listOfPairs.size() || index.row() < 0) + return QVariant(); + + if (role == Qt::DisplayRole) { + QPair pair = listOfPairs.at(index.row()); + + if (index.column() == 0) + return pair.first; + else if (index.column() == 1) + return pair.second; + } + return QVariant(); +} + +QVariant QmitkTbssMetaTableModel::headerData(int section, Qt::Orientation orientation, int role) const +{ + if (role != Qt::DisplayRole) + return QVariant(); + + if (orientation == Qt::Horizontal) { + switch (section) { + case 0: + return tr("Image"); + + case 1: + return tr("Filename"); + + default: + return QVariant(); + } + } + return QVariant(); +} + +bool QmitkTbssMetaTableModel::insertRows(int position, int rows, const QModelIndex &index) +{ + Q_UNUSED(index); + beginInsertRows(QModelIndex(), position, position+rows-1); + + for (int row=0; row < rows; row++) { + QPair pair("Image ", "Filename"); + listOfPairs.insert(position, pair); + } + + endInsertRows(); + return true; +} + +bool QmitkTbssMetaTableModel::removeRows(int position, int rows, const QModelIndex &index) +{ + Q_UNUSED(index); + beginRemoveRows(QModelIndex(), position, position+rows-1); + + for (int row=0; row < rows; ++row) { + listOfPairs.removeAt(position); + } + + endRemoveRows(); + return true; +} + +bool QmitkTbssMetaTableModel::setData(const QModelIndex &index, const QVariant &value, int role) +{ + if (index.isValid() && role == Qt::EditRole) { + int row = index.row(); + + QPair p = listOfPairs.value(row); + + if (index.column() == 0) + p.first = value.toString(); + else if (index.column() == 1) + p.second = value.toString(); + else + return false; + + listOfPairs.replace(row, p); + emit(dataChanged(index, index)); + + return true; + } + + return false; +} + +Qt::ItemFlags QmitkTbssMetaTableModel::flags(const QModelIndex &index) const +{ + if (!index.isValid()) + return Qt::ItemIsEnabled; + + return QAbstractTableModel::flags(index) | Qt::ItemIsEditable; +} + +QList< QPair > QmitkTbssMetaTableModel::getList() +{ + return listOfPairs; +} + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.h new file mode 100644 index 0000000000..112b33bc45 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.h @@ -0,0 +1,55 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +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 QmitkTbssMetaTableModel_h +#define QmitkTbssMetaTableModel_h + + +#include + +#include +#include +#include +#include + + +class QmitkTbssMetaTableModel : public QAbstractTableModel +{ + + //Q_OBJECT + +public: + QmitkTbssMetaTableModel(QObject *parent=0); + QmitkTbssMetaTableModel(QList< QPair > listofPairs, QObject *parent=0); + + int rowCount(const QModelIndex &parent) const; + int columnCount(const QModelIndex &parent) const; + QVariant data(const QModelIndex &index, int role) const; + QVariant headerData(int section, Qt::Orientation orientation, int role) const; + Qt::ItemFlags flags(const QModelIndex &index) const; + bool setData(const QModelIndex &index, const QVariant &value, int role=Qt::EditRole); + bool insertRows(int position, int rows, const QModelIndex &index=QModelIndex()); + bool removeRows(int position, int rows, const QModelIndex &index=QModelIndex()); + QList< QPair > getList(); + +private: + QList< QPair > listOfPairs; +}; + +#endif diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssTableModel.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssTableModel.cpp new file mode 100644 index 0000000000..6319490fff --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssTableModel.cpp @@ -0,0 +1,150 @@ + + +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +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. + +=========================================================================*/ + + + +#include "QmitkTbssTableModel.h" + +QmitkTbssTableModel::QmitkTbssTableModel(QObject *parent) + : QAbstractTableModel(parent) +{ +} + +QmitkTbssTableModel::QmitkTbssTableModel(QList< QPair > pairs, QObject *parent) + : QAbstractTableModel(parent) +{ + listOfPairs=pairs; +} + +int QmitkTbssTableModel::rowCount(const QModelIndex &parent) const +{ + Q_UNUSED(parent); + return listOfPairs.size(); +} + +int QmitkTbssTableModel::columnCount(const QModelIndex &parent) const +{ + Q_UNUSED(parent); + return 2; +} + +QVariant QmitkTbssTableModel::data(const QModelIndex &index, int role) const +{ + if (!index.isValid()) + return QVariant(); + + if (index.row() >= listOfPairs.size() || index.row() < 0) + return QVariant(); + + if (role == Qt::DisplayRole) { + QPair pair = listOfPairs.at(index.row()); + + if (index.column() == 0) + return pair.first; + else if (index.column() == 1) + return pair.second; + } + return QVariant(); +} + +QVariant QmitkTbssTableModel::headerData(int section, Qt::Orientation orientation, int role) const +{ + if (role != Qt::DisplayRole) + return QVariant(); + + if (orientation == Qt::Horizontal) { + switch (section) { + case 0: + return tr("Group"); + + case 1: + return tr("Size"); + + default: + return QVariant(); + } + } + return QVariant(); +} + +bool QmitkTbssTableModel::insertRows(int position, int rows, const QModelIndex &index) +{ + Q_UNUSED(index); + beginInsertRows(QModelIndex(), position, position+rows-1); + + for (int row=0; row < rows; row++) { + QPair pair(" ", 0); + listOfPairs.insert(position, pair); + } + + endInsertRows(); + return true; +} + +bool QmitkTbssTableModel::removeRows(int position, int rows, const QModelIndex &index) +{ + Q_UNUSED(index); + beginRemoveRows(QModelIndex(), position, position+rows-1); + + for (int row=0; row < rows; ++row) { + listOfPairs.removeAt(position); + } + + endRemoveRows(); + return true; +} + +bool QmitkTbssTableModel::setData(const QModelIndex &index, const QVariant &value, int role) +{ + if (index.isValid() && role == Qt::EditRole) { + int row = index.row(); + + QPair p = listOfPairs.value(row); + + if (index.column() == 0) + p.first = value.toString(); + else if (index.column() == 1) + p.second = value.toInt(); + else + return false; + + listOfPairs.replace(row, p); + emit(dataChanged(index, index)); + + return true; + } + + return false; +} + +Qt::ItemFlags QmitkTbssTableModel::flags(const QModelIndex &index) const +{ + if (!index.isValid()) + return Qt::ItemIsEnabled; + + return QAbstractTableModel::flags(index) | Qt::ItemIsEditable; +} + +QList< QPair > QmitkTbssTableModel::getList() +{ + return listOfPairs; +} + + + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssTableModel.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssTableModel.h new file mode 100644 index 0000000000..3f900fd02b --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssTableModel.h @@ -0,0 +1,53 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date$ +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 QmitkTbssTableModel_h +#define QmitkTbssTableModel_h + + +#include + +#include +#include +#include + +class QmitkTbssTableModel : public QAbstractTableModel +{ + + //Q_OBJECT + +public: + QmitkTbssTableModel(QObject *parent=0); + QmitkTbssTableModel(QList< QPair > listofPairs, QObject *parent=0); + + int rowCount(const QModelIndex &parent) const; + int columnCount(const QModelIndex &parent) const; + QVariant data(const QModelIndex &index, int role) const; + QVariant headerData(int section, Qt::Orientation orientation, int role) const; + Qt::ItemFlags flags(const QModelIndex &index) const; + bool setData(const QModelIndex &index, const QVariant &value, int role=Qt::EditRole); + bool insertRows(int position, int rows, const QModelIndex &index=QModelIndex()); + bool removeRows(int position, int rows, const QModelIndex &index=QModelIndex()); + QList< QPair > getList(); + +private: + QList< QPair > listOfPairs; +}; + +#endif diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp new file mode 100644 index 0000000000..0a5be3ac4a --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp @@ -0,0 +1,1601 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ +Version: $Revision: 21975 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +// Blueberry +#include "berryIWorkbenchWindow.h" +#include "berryIWorkbenchPage.h" +#include "berryISelectionService.h" +#include "berryConstants.h" +#include "berryPlatformUI.h" + +// Qmitk +#include "QmitkTractbasedSpatialStatisticsView.h" +#include "QmitkStdMultiWidget.h" + +#include "mitkDataNodeObject.h" +#include + +// Qt +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +//#include +//#include +//#include +//#include +#include +//#include +//#include +//#include +//#include +//#include +//#include +#include +//#include +#include + + +#include "mitkPartialVolumeAnalysisClusteringCalculator.h" +#include "mitkPartialVolumeAnalysisHistogramCalculator.h" + +#include +#include "vtkFloatArray.h" +#include "vtkLinearTransform.h" +#include "vtkPoints.h" +#include "mitkSurface.h" +#include +#include "vtkArrowSource.h" +#include "vtkUnstructuredGrid.h" +#include "vtkPointData.h" + +#include +#include +#include + +#include + +#include "mitkITKImageImport.h" +// #include "mitkImageMapperGL2D.h" +#include "mitkVolumeDataVtkMapper3D.h" +#include "mitkImageAccessByItk.h" + + +#define SEARCHSIGMA 10 /* length in linear voxel dimens + { + // create new ones + m_PointSetNode = mitk::PointSet::New();ions */ +#define MAXSEARCHLENGTH (3*SEARCHSIGMA) + +const std::string QmitkTractbasedSpatialStatisticsView::VIEW_ID = "org.mitk.views.tractbasedspatialstatistics"; + +using namespace berry; + + +struct TbssSelListener : ISelectionListener +{ + + berryObjectMacro(TbssSelListener) + + TbssSelListener(QmitkTractbasedSpatialStatisticsView* view) + { + m_View = view; + } + + + void DoSelectionChanged(ISelection::ConstPointer selection) + { + // save current selection in member variable + m_View->m_CurrentSelection = selection.Cast(); + + // do something with the selected items + if(m_View->m_CurrentSelection) + { + bool foundTbssRoi = false; + bool foundTbss = false; + bool found3dImage = false; + bool found4dImage = false; + + mitk::TbssRoiImage* roiImage; + mitk::TbssImage* image; + + + // iterate selection + for (IStructuredSelection::iterator i = m_View->m_CurrentSelection->Begin(); + i != m_View->m_CurrentSelection->End(); ++i) + { + + // extract datatree node + if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) + { + mitk::DataNode::Pointer node = nodeObj->GetDataNode(); + + // only look at interesting types + + if(QString("TbssRoiImage").compare(node->GetData()->GetNameOfClass())==0) + { + foundTbssRoi = true; + roiImage = static_cast(node->GetData()); + } + else if (QString("TbssImage").compare(node->GetData()->GetNameOfClass())==0) + { + foundTbss = true; + image = static_cast(node->GetData()); + } + else if(QString("Image").compare(node->GetData()->GetNameOfClass())==0) + { + mitk::Image* img = static_cast(node->GetData()); + if(img->GetDimension() == 3) + { + found3dImage = true; + } + else if(img->GetDimension() == 4) + { + found4dImage = true; + } + } + + + } + + } + + + m_View->m_Controls->m_CreateRoi->setEnabled(found3dImage); + m_View->m_Controls->m_ImportFsl->setEnabled(found4dImage); + if(found3dImage) + { + m_View->InitPointsets(); + } + + } + } + + + void SelectionChanged(IWorkbenchPart::Pointer part, ISelection::ConstPointer selection) + { + // check, if selection comes from datamanager + if (part) + { + QString partname(part->GetPartName().c_str()); + if(partname.compare("Datamanager")==0) + { + // apply selection + DoSelectionChanged(selection); + } + } + } + + QmitkTractbasedSpatialStatisticsView* m_View; +}; + + +QmitkTractbasedSpatialStatisticsView::QmitkTractbasedSpatialStatisticsView() +: QmitkFunctionality() +, m_Controls( 0 ) +, m_MultiWidget( NULL ) +{ + +} + +QmitkTractbasedSpatialStatisticsView::~QmitkTractbasedSpatialStatisticsView() +{ +} + +void QmitkTractbasedSpatialStatisticsView::OnSelectionChanged(std::vector nodes) +{ + //datamanager selection changed + if (!this->IsActivated()) + return; + + + // Get DataManagerSelection + if (!this->GetDataManagerSelection().empty()) + { + mitk::DataNode::Pointer sourceImageNode = this->GetDataManagerSelection().front(); + mitk::Image::Pointer sourceImage = dynamic_cast(sourceImageNode->GetData()); + + if (!sourceImage) + { + m_Controls->m_TbssImageLabel->setText( + QString( sourceImageNode->GetName().c_str() ) + " is no image" + ); + + return; + } + + // set Text + m_Controls->m_TbssImageLabel->setText( + QString( sourceImageNode->GetName().c_str() ) + " (" + + QString::number(sourceImage->GetDimension()) + "D)" + ); + + + } + else + { + m_Controls->m_TbssImageLabel->setText("Please select an image"); + } + +} + +void QmitkTractbasedSpatialStatisticsView::InitPointsets() +{ + // Check if PointSetStart exsits, if not create it. + m_P1 = this->GetDefaultDataStorage()->GetNamedNode("PointSetNode"); + + if (m_PointSetNode) + { + //m_PointSetNode = dynamic_cast(m_P1->GetData()); + return; + } + + if ((!m_P1) || (!m_PointSetNode)) + { + // create new ones + m_PointSetNode = mitk::PointSet::New(); + m_P1 = mitk::DataNode::New(); + m_P1->SetData( m_PointSetNode ); + m_P1->SetProperty( "name", mitk::StringProperty::New( "PointSet" ) ); + m_P1->SetProperty( "opacity", mitk::FloatProperty::New( 1 ) ); + m_P1->SetProperty( "helper object", mitk::BoolProperty::New(false) ); // CHANGE if wanted + m_P1->SetProperty( "pointsize", mitk::FloatProperty::New( 0.1 ) ); + m_P1->SetColor( 1.0, 0.0, 0.0 ); + this->GetDefaultDataStorage()->Add(m_P1); + m_Controls->m_PointWidget->SetPointSetNode(m_P1); + m_Controls->m_PointWidget->SetMultiWidget(GetActiveStdMultiWidget()); + } +} + +void QmitkTractbasedSpatialStatisticsView::CreateQtPartControl( QWidget *parent ) +{ + // build up qt view, unless already done + if ( !m_Controls ) + { + // create GUI widgets from the Qt Designer's .ui file + m_Controls = new Ui::QmitkTractbasedSpatialStatisticsViewControls; + m_Controls->setupUi( parent ); + this->CreateConnections(); + } + + m_SelListener = berry::ISelectionListener::Pointer(new TbssSelListener(this)); + this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->AddPostSelectionListener(/*"org.mitk.views.datamanager",*/ m_SelListener); + berry::ISelection::ConstPointer sel( + this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); + m_CurrentSelection = sel.Cast(); + m_SelListener.Cast()->DoSelectionChanged(sel); + + m_IsInitialized = false; + + + + // Table for the FSL TBSS import + m_GroupModel = new QmitkTbssTableModel(); + m_Controls->m_GroupInfo->setModel(m_GroupModel); + + +} + +void QmitkTractbasedSpatialStatisticsView::Activated() +{ + QmitkFunctionality::Activated(); + + berry::ISelection::ConstPointer sel( + this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); + m_CurrentSelection = sel.Cast(); + m_SelListener.Cast()->DoSelectionChanged(sel); +} + +void QmitkTractbasedSpatialStatisticsView::Deactivated() +{ + QmitkFunctionality::Deactivated(); +} + +void QmitkTractbasedSpatialStatisticsView::CreateConnections() +{ + if ( m_Controls ) + { + connect( (QObject*)(m_Controls->m_CreateRoi), SIGNAL(clicked()), this, SLOT(CreateRoi()) ); + connect( (QObject*)(m_Controls->m_ImportFsl), SIGNAL(clicked()), this, SLOT(TbssImport()) ); + connect( (QObject*)(m_Controls->m_AddGroup), SIGNAL(clicked()), this, SLOT(AddGroup()) ); + connect( (QObject*)(m_Controls->m_RemoveGroup), SIGNAL(clicked()), this, SLOT(RemoveGroup()) ); + connect( (QObject*)(m_Controls->m_Clipboard), SIGNAL(clicked()), this, SLOT(CopyToClipboard()) ); + connect( m_Controls->m_RoiPlotWidget->m_PlotPicker, SIGNAL(selected(const QwtDoublePoint&)), SLOT(Clicked(const QwtDoublePoint&) ) ); + connect( m_Controls->m_RoiPlotWidget->m_PlotPicker, SIGNAL(moved(const QwtDoublePoint&)), SLOT(Clicked(const QwtDoublePoint&) ) ); + } +} + +void QmitkTractbasedSpatialStatisticsView::CopyToClipboard() +{ + std::vector > vals = m_Controls->m_RoiPlotWidget->GetVals(); + 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)); + + double d = *it2; + std::cout << d <setText(clipboardText, QClipboard::Clipboard); + +} + +void QmitkTractbasedSpatialStatisticsView::RemoveGroup() +{ + + QTableView *temp = static_cast(m_Controls->m_GroupInfo); + // QSortFilterProxyModel *proxy = static_cast(temp->model()); + QItemSelectionModel *selectionModel = temp->selectionModel(); + + QModelIndexList indices = selectionModel->selectedRows(); + + QModelIndex index; + + foreach(index, indices) + { + int row = index.row(); + m_GroupModel->removeRows(row, 1, QModelIndex()); + } + +} + +std::string QmitkTractbasedSpatialStatisticsView::ReadFile(std::string whatfile) +{ + std::string s = "Select a" + whatfile; + QFileDialog* w = new QFileDialog(this->m_Controls->m_ImportFsl, QString(s.c_str()) ); + w->setFileMode(QFileDialog::ExistingFiles); + w->setDirectory("/home"); + + if(whatfile == "gradient image") + { + w->setNameFilter("Tbss gradient images (*.tgi)"); + } + + // RETRIEVE SELECTION + if ( w->exec() != QDialog::Accepted ) + { + return ""; + MITK_INFO << "Failed to load"; + } + + QStringList filenames = w->selectedFiles(); + if (filenames.size() > 0) + { + std::string retval = filenames.at(0).toStdString(); + return retval; + } + return ""; +} + +void QmitkTractbasedSpatialStatisticsView::AddGroup() +{ + QString group("Group"); + int number = 0; + QPair pair(group, number); + QList< QPair >list = m_GroupModel->getList(); + + + if(!list.contains(pair)) + { + m_GroupModel->insertRows(0, 1, QModelIndex()); + + QModelIndex index = m_GroupModel->index(0, 0, QModelIndex()); + m_GroupModel->setData(index, group, Qt::EditRole); + index = m_GroupModel->index(0, 1, QModelIndex()); + m_GroupModel->setData(index, number, Qt::EditRole); + + } + else + { + //QMessageBox::information(this, "Duplicate name"); + } + + +} + + +void QmitkTractbasedSpatialStatisticsView::TbssImport() +{ + + // Read groups from the interface + mitk::TbssImporter::Pointer importer = mitk::TbssImporter::New(); + + QList< QPair >list = m_GroupModel->getList(); + + if(list.size() == 0) + { + QMessageBox msgBox; + msgBox.setText("No study group information has been set yet."); + msgBox.exec(); + return; + } + + + std::vector < std::pair > groups; + for(int i=0; i pair = list.at(i); + std::string s = pair.first.toStdString(); + int n = pair.second; + + std::pair p; + p.first = s; + p.second = n; + groups.push_back(p); + } + + importer->SetGroupInfo(groups); + + std::string minfo = m_Controls->m_MeasurementInfo->text().toStdString(); + importer->SetMeasurementInfo(minfo); + + + std::string name = ""; + for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); + i != m_CurrentSelection->End(); ++i) + { + // extract datatree node + if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) + { + mitk::DataNode::Pointer node = nodeObj->GetDataNode(); + + if(QString("Image").compare(node->GetData()->GetNameOfClass())==0) + { + mitk::Image* img = static_cast(node->GetData()); + if(img->GetDimension() == 4) + { + //importer->SetImportVolume(img); + name = node->GetName(); + } + } + } + } + + mitk::TbssImage::Pointer tbssImage; + + tbssImage = importer->Import(); + name += "_tbss"; + AddTbssToDataStorage(tbssImage, name); + + +} + + +void QmitkTractbasedSpatialStatisticsView::AddTbssToDataStorage(mitk::Image* image, std::string name) +{ + mitk::LevelWindow levelwindow; + levelwindow.SetAuto( image ); + mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); + levWinProp->SetLevelWindow( levelwindow ); + + mitk::DataNode::Pointer result = mitk::DataNode::New(); + result->SetProperty( "name", mitk::StringProperty::New(name) ); + result->SetData( image ); + result->SetProperty( "levelwindow", levWinProp ); + + + // add new image to data storage and set as active to ease further processing + GetDefaultDataStorage()->Add( result ); + + // show the results + mitk::RenderingManager::GetInstance()->RequestUpdateAll(); +} + +void QmitkTractbasedSpatialStatisticsView::Clicked(const QwtDoublePoint& pos) +{ + if(m_Roi.size() > 0 && m_CurrentGeometry != NULL) + { + int index = (int)pos.x(); + index = std::max(0, index); + index = std::min(index, (int)m_Roi.size()); + itk::Index<3> ix = m_Roi.at(index); + + mitk::Vector3D i; + i[0] = ix[0]; + i[1] = ix[1]; + i[2] = ix[2]; + + mitk::Vector3D w; + + m_CurrentGeometry->IndexToWorld(i, w); + + mitk::Point3D origin = m_CurrentGeometry->GetOrigin(); + + mitk::Point3D p; + p[0] = w[0] + origin[0]; + p[1] = w[1] + origin[1]; + p[2] = w[2] + origin[2]; + + + + + m_MultiWidget->MoveCrossToPosition(p); + + m_Controls->m_RoiPlotWidget->drawBar(index); + + } + +} + +void QmitkTractbasedSpatialStatisticsView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) +{ + m_MultiWidget = &stdMultiWidget; +} + + +void QmitkTractbasedSpatialStatisticsView::StdMultiWidgetNotAvailable() +{ + m_MultiWidget = NULL; +} + +void QmitkTractbasedSpatialStatisticsView::AdjustPlotMeasure(const QString & text) +{ + berry::ISelection::ConstPointer sel( + this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); + m_CurrentSelection = sel.Cast(); + m_SelListener.Cast()->DoSelectionChanged(sel); +} + +void QmitkTractbasedSpatialStatisticsView::Clustering() +{ +/* + + + // Create a mask using the distance map + typedef itk::ImageFileReader< VectorImageType > DirectionReader; + DirectionReader::Pointer directionReader = DirectionReader::New(); + directionReader->SetFileName(m_TbssWorkspaceManager.GetInputDir().toStdString() + "/tbss/" + m_TbssWorkspaceManager.GetGradient().toStdString()); + directionReader->Update(); + VectorImageType::Pointer directions = directionReader->GetOutput(); + + FloatReaderType::Pointer distMapReader = FloatReaderType::New(); + distMapReader->SetFileName(m_TbssWorkspaceManager.GetInputDir().toStdString() + "/stats/" + m_TbssWorkspaceManager.GetDistanceMap().toStdString()); + distMapReader->Update(); + FloatImageType::Pointer distanceMap = distMapReader->GetOutput(); + + + std::string line; + std::string path = "/mnt/E130-Projekte/NeuroDiffusion/BRAM DTI/TBSS/rois/cc.txt"; + std::ifstream file(path.c_str()); + std::vector< itk::Index< 3 > > roi; + + if(file.is_open()) + { + + while(getline(file,line)) + { + + std::vector tokens; + Tokenize(line, tokens); + + itk::Index<3> ix; + ix[0] = atoi(tokens[0].c_str()); + ix[1] = atoi(tokens[1].c_str()); + ix[2] = atoi(tokens[2].c_str()); + roi.push_back(ix); + } + } + + if(roi.size() == 0) + { + return; + } + + + + // Some code from the projection algorithm of tbss to create a mask + + + std::vector< std::vector< itk::Index< 3 > > > rois; + + for(int j=0; j > indices; + + FloatImageType::SizeType size = distanceMap->GetLargestPossibleRegion().GetSize(); + + bool roiDone = false; + + while(!roiDone && j ix = roi[j]; + int x=ix[0]; + int y=ix[1]; + int z=ix[2]; + VectorImageType::PixelType dir = directions->GetPixel(ix); + + indices.push_back(ix); + + + for(int iters=0;iters<2;iters++) + { + float distance=0; + + for(int d=1;d=size[0] && dy<=size[1] && dz<=size[2]) + { + d=MAXSEARCHLENGTH; + } + else if(distanceMap->GetPixel(ix)>=distance) + { + distance = distanceMap->GetPixel(ix); + indices.push_back(ix); + } + else{ + d=MAXSEARCHLENGTH; + } + + } + + } + j++; + } + + + // Create a mask from indices + UCharImageType::Pointer maskItk = UCharImageType::New(); + maskItk->SetRegions(distanceMap->GetRequestedRegion()); + maskItk->SetDirection(distanceMap->GetDirection()); + maskItk->SetSpacing(distanceMap->GetSpacing()); + maskItk->SetOrigin(distanceMap->GetOrigin()); + maskItk->Allocate(); + + + + + + + + // For every point on the roi create a mask and feed it to the partial voluming algorithm + + + //maskItk->FillBuffer(0); + + + + // Create a bounding box from current ROI + int xMin = numeric_limits::max(); + int yMin = numeric_limits::max(); + int zMin = numeric_limits::max(); + int xMax = numeric_limits::min(); + int yMax = numeric_limits::min(); + int zMax = numeriUCharImageType::Pointer newMask = UCharImageType::New();c_limits::min(); + + + for(int i=0; i ix = indices[i]; + + if(ix[0] < xMin) + xMin=ix[0]; + + if(ix[1] < yMin) + yMin=ix[1]; + + if(ix[2] < zMin) + zMin=ix[2]; + + if(ix[0] > xMax) + xMax=ix[0]; + + if(ix[1] > yMax) + yMax=ix[1]; + + if(ix[2] > zMax) + zMax=ix[2]; + } + + FloatImageType::PointType origin = distanceMap->GetOrigin(); + CharImageType::PointType originMask; + originMask[0] = origin[0] + xMin; + originMask[1] = origin[1] + -yMin; + originMask[2] = origin[2] + zMin; + + CharImageType::RegionType region; + CharImageType::RegionType::SizeType s; + s[0] = xMax-xMin + 1; + s[1] = yMax-yMin + 1; + s[2] = zMax-zMin + 1; + region.SetSize(s); + + + UCharImageType::Pointer newMask = UCharImageType::New(); + newMask->SetSpacing( distanceMap->GetSpacing() ); // Set the image spacing + newMask->SetOrigin( originMask ); // Set the image origin + newMask->SetDirection( distanceMap->GetDirection() ); // Set the image direction + newMask->SetRegions( region ); + newMask->Allocate(); + newMask->FillBuffer(0); + + + + for(int i=0; i ix = indices[i]; + + itk::Point< double, 3 > point; + itk::Index< 3 > index; + + distanceMap->TransformIndexToPhysicalPoint (ix, point); + + + + newMask->TransformPhysicalPointToIndex(point, index); + + + + newMask->SetPixel(index, 1); + } + + + +*/ + + + + + + + + UCharImageType::Pointer newMask = UCharImageType::New(); + UCharReaderType::Pointer cReader = UCharReaderType::New(); + cReader->SetFileName("/mnt/E130-Projekte/NeuroDiffusion/BRAM DTI/ClusteringFornix/fornix_central_maxFA_path_Dilated_by_3.nrrd"); + cReader->Update(); + newMask = cReader->GetOutput(); + + + // mitk::DataNode::Pointer maskNode = readNode("itk image/mnt/E130-Projekte/NeuroDiffusion/BRAM DTI/TBSS/clusterMasks/area2.nii"); + // mitk::Image::Pointer mask = dynamic_cast(maskNode->GetData()); + mitk::Image::Pointer mask; + mitk::CastToMitkImage(newMask, mask); + + typedef mitk::PartialVolumeAnalysisHistogramCalculator HistorgramCalculator; + typedef mitk::PartialVolumeAnalysisClusteringCalculator ClusteringType; + typedef HistorgramCalculator::HistogramType HistogramType; + + HistorgramCalculator::Pointer histogramCalculator = HistorgramCalculator::New(); + + + + + + // Make list of subjects + + std::vector paths; + paths.push_back("/mnt/E130-Projekte/NeuroDiffusion/BRAM DTI/TBSS/FA/SORTEDBYCONDITION/FA/subset"); + // paths.push_back("/mnt/E130-Projekte/NeuroDiffusion/BRAM DTI/TBSS/FA/SORTEDBYCONDITION/AXD/"); + + for(int j=0; j values; + + for(int i=0; i(node->GetData()); + + + histogramCalculator->SetImage(image); + histogramCalculator->SetImageMask( mask ); + histogramCalculator->SetMaskingModeToImage(); + histogramCalculator->SetNumberOfBins(25); + histogramCalculator->SetUpsamplingFactor(5); + histogramCalculator->SetGaussianSigma(0.0); + histogramCalculator->SetForceUpdate(true); + + bool statisticsChanged = histogramCalculator->ComputeStatistics( ); + + + ClusteringType::ParamsType *cparams = 0; + ClusteringType::ClusterResultType *cresult = 0; + ClusteringType::HistType *chist = 0; + ClusteringType::HelperStructPerformClusteringRetval *currentPerformClusteringResults; + + try{ + + + mitk::Image* tmpImg = histogramCalculator->GetInternalImage(); + mitk::Image::ConstPointer imgToCluster = tmpImg; + + if(imgToCluster.IsNotNull()) + { + + // perform clustering + const HistogramType *histogram = histogramCalculator->GetHistogram( ); + ClusteringType::Pointer clusterer = ClusteringType::New(); + clusterer->SetStepsNumIntegration(200); + clusterer->SetMaxIt(1000); + + mitk::Image::Pointer pFiberImg; + + currentPerformClusteringResults = + clusterer->PerformClustering(imgToCluster, histogram, 2); + + pFiberImg = currentPerformClusteringResults->clusteredImage; + cparams = currentPerformClusteringResults->params; + cresult = currentPerformClusteringResults->result; + chist = currentPerformClusteringResults->hist; + + + // m_Controls->m_HistogramWidget->SetParameters( + // cparams, cresult, chist ); + + + + std::vector *xVals = chist->GetXVals(); + + std::vector *fiberVals = new std::vector(cresult->GetFiberVals()); + + + + double fiberFA = 0.0; + double weights = 0.0; + + // std::cout << "x, y, fiber, nonFiber, mixed, combi" << std::endl; + for(int k=0; ksize(); ++k) + { + + fiberFA += xVals->at(k) * fiberVals->at(k); + weights += fiberVals->at(k); + + + } + + fiberFA = fiberFA / weights; + std::cout << "FA: " << fiberFA << std::endl; + values.push_back(fiberFA); + + } + + + } + + catch ( const std::runtime_error &e ) + { + std::cout << "noooooooooooooooooooooooooooooooo!"; + } + + //MITK_INFO << "number of voxels: " << indices.size(); + } + + std::vector::iterator it = values.begin(); + while(it!=values.end()) + { + std::cout << *it << std::endl; + ++it; + } + + } +} + + +void QmitkTractbasedSpatialStatisticsView::CreateRoi() +{ + + // It is important to load the MeanFASkeletonMask image in MITK to make sure that point selection and + // pathfinding is done on the same image + //string filename = m_TbssWorkspaceManager.GetInputDir().toStdString() + "/stats/" + m_TbssWorkspaceManager.GetMeanFASkeletonMask().toStdString(); + + // Implement a way to obtain skeleton and skeletonFA without sml workspace + + double threshold = QInputDialog::getDouble(m_Controls->m_CreateRoi, tr("Set an FA threshold"), + tr("Threshold:"), QLineEdit::Normal, + 0.2); + + + mitk::Image::Pointer image; + + for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); + i != m_CurrentSelection->End(); ++i) + { + // extract datatree node + if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) + { + mitk::DataNode::Pointer node = nodeObj->GetDataNode(); + + if(QString("Image").compare(node->GetData()->GetNameOfClass())==0) + { + mitk::Image* img = static_cast(node->GetData()); + if(img->GetDimension() == 3) + { + image = img; + } + } + } + } + + + + if(image.IsNull()) + { + return; + } + + mitk::TractAnalyzer analyzer; + analyzer.SetInputImage(image); + analyzer.SetThreshold(threshold); + + + + int n = 0; + if(m_PointSetNode.IsNotNull()) + { + n = m_PointSetNode->GetSize(); + if(n==0) + { + QMessageBox msgBox; + msgBox.setText("No points have been set yet."); + msgBox.exec(); + } + } + else{ + QMessageBox msgBox; + msgBox.setText("No points have been set yet."); + msgBox.exec(); + } + + std::string pathDescription = ""; + std::vector< itk::Index<3> > totalPath; + + if(n>0) + { + for(int i=0; iGetPoint(i); + mitk::Point3D p2 = m_PointSetNode->GetPoint(i+1); + + + itk::Index<3> StartPoint; + itk::Index<3> EndPoint; + image->GetGeometry()->WorldToIndex(p,StartPoint); + image->GetGeometry()->WorldToIndex(p2,EndPoint); + MITK_INFO << "create roi"; + + + analyzer.BuildGraph(StartPoint, EndPoint); + std::vector< itk::Index<3> > path = analyzer.GetPath(); + + + for(std::vector< itk::Index<3> >::iterator it = path.begin(); + it != path.end(); it++) + { + itk::Index<3> ix = *it; + + if (!(ix==EndPoint)) + { + totalPath.push_back(ix); + std::stringstream ss; + ss << ix[0] << " " << ix[1] << " " << ix[2] << "\n"; + pathDescription += ss.str(); + } + else + { + // Only when dealing with the last segment the last point should be added. This one will not occur + // as the first point of the next roi segment. + if(i == (n-2)) + { + totalPath.push_back(EndPoint); + std::stringstream ss; + ss << EndPoint[0] << " " << EndPoint[1] << " " << EndPoint[2] << "\n"; + pathDescription += ss.str(); + } + + } + } + + + + } + + m_Controls->m_PathTextEdit->setPlainText(QString(pathDescription.c_str())); + + + FloatImageType::Pointer itkImg = FloatImageType::New(); + mitk::CastToItkImage(image, itkImg); + + CharImageType::Pointer roiImg = CharImageType::New(); + roiImg->SetRegions(itkImg->GetLargestPossibleRegion().GetSize()); + roiImg->SetOrigin(itkImg->GetOrigin()); + roiImg->SetSpacing(itkImg->GetSpacing()); + roiImg->SetDirection(itkImg->GetDirection()); + roiImg->Allocate(); + roiImg->FillBuffer(0); + + + std::vector< itk::Index<3> > roi; + + std::vector< itk::Index<3> >::iterator it; + for(it = totalPath.begin(); + it != totalPath.end(); + it++) + { + itk::Index<3> ix = *it; + roiImg->SetPixel(ix, 1); + roi.push_back(ix); + } + + + mitk::TbssRoiImage::Pointer tbssRoi = mitk::TbssRoiImage::New(); + //mitk::CastToTbssImage(m_CurrentRoi.GetPointer(), tbssRoi); + tbssRoi->SetRoi(roi); + tbssRoi->SetImage(roiImg); + tbssRoi->SetStructure(m_Controls->m_Structure->text().toStdString()); + tbssRoi->InitializeFromImage(); + + + // mitk::Image::Pointer tbssRoi = mitk::Image::New(); + //mitk::CastToTbssImage(m_CurrentRoi.GetPointer(), tbssRoi); + // mitk::CastToMitkImage(roiImg, tbssRoi); + + AddTbssToDataStorage(tbssRoi, m_Controls->m_RoiName->text().toStdString()); + + } + +} + + + + + +void QmitkTractbasedSpatialStatisticsView::Plot(mitk::TbssImage* image, mitk::TbssRoiImage* roiImage) +{ + if(m_Controls->m_TabWidget->currentWidget() == m_Controls->m_MeasureTAB) + { + + std::vector< itk::Index<3> > roi = roiImage->GetRoi(); + m_Roi = roi; + m_CurrentGeometry = image->GetGeometry(); + + + std::string resultfile = ""; + + /* + if(image->GetPreprocessedFA()) + { + resultFile = image->GetPreprocessedFAFile(); + } + */ + std::string structure = roiImage->GetStructure(); + + + + //m_View->m_CurrentGeometry = image->GetGeometry(); + + m_Controls->m_RoiPlotWidget->SetGroups(image->GetGroupInfo()); + + + // Check for preprocessed results to save time + + //if(resultfile == "") + // { + // Need to calculate the results using the 4D volume + // Can save the time this takes if there are results available already + + //std::string type = m_Controls->m_MeasureType->itemText(m_Controls->m_MeasureType->currentIndex()).toStdString(); + m_Controls->m_RoiPlotWidget->SetProjections(image->GetImage()); + + + // } + + m_Controls->m_RoiPlotWidget->SetRoi(roi); + m_Controls->m_RoiPlotWidget->SetStructure(structure); + m_Controls->m_RoiPlotWidget->SetMeasure( image->GetMeasurementInfo() ); + m_Controls->m_RoiPlotWidget->DrawProfiles(resultfile); + } +} + + + + +void QmitkTractbasedSpatialStatisticsView::Masking() +{ + //QString filename = m_Controls->m_WorkingDirectory->text(); + QString filename = "E:/Experiments/tbss"; + QString faFiles = filename + "/AxD"; + QString maskFiles = filename + "/bin_masks"; + + + QDirIterator faDirIt(faFiles, QDir::Files | QDir::NoSymLinks, QDirIterator::Subdirectories); + QDirIterator maskDirIt(maskFiles, QDir::Files | QDir::NoSymLinks, QDirIterator::Subdirectories); + + std::vector faFilenames; + std::vector maskFilenames; + std::vector outputFilenames; + + while(faDirIt.hasNext() && maskDirIt.hasNext()) + { + faDirIt.next(); + maskDirIt.next(); + if((faDirIt.fileInfo().completeSuffix() == "nii" + || faDirIt.fileInfo().completeSuffix() == "mhd" + || faDirIt.fileInfo().completeSuffix() == "nii.gz") + && (maskDirIt.fileInfo().completeSuffix() == "nii" + || maskDirIt.fileInfo().completeSuffix() == "mhd" + || maskDirIt.fileInfo().completeSuffix() == "nii.gz")) + { + faFilenames.push_back(faDirIt.filePath().toStdString()); + outputFilenames.push_back(faDirIt.fileName().toStdString()); + maskFilenames.push_back(maskDirIt.filePath().toStdString()); + } + } + + std::vector::iterator faIt = faFilenames.begin(); + std::vector::iterator maskIt = maskFilenames.begin(); + std::vector::iterator outputIt = outputFilenames.begin(); + + // Now multiply all FA images with their corresponding masks + + QString outputDir = filename; + while(faIt != faFilenames.end() && maskIt != maskFilenames.end() && outputIt != outputFilenames.end()) + { + std::cout << "Mask " << *faIt << " with " << *maskIt << std::endl; + + typedef itk::MultiplyImageFilter MultiplicationFilterType; + + FloatReaderType::Pointer floatReader = FloatReaderType::New(); + CharReaderType::Pointer charReader = CharReaderType::New(); + + floatReader->SetFileName(*faIt); + //floatReader->Update(); + //FloatImageType::Pointer faImage = floatReader->GetOutput(); + + charReader->SetFileName(*maskIt); + //charReader->Update(); + // CharImageType::Pointer maskImage = charReader->GetOutput(); + + MultiplicationFilterType::Pointer multiplicationFilter = MultiplicationFilterType::New(); + multiplicationFilter->SetInput1(floatReader->GetOutput()); + multiplicationFilter->SetInput2(charReader->GetOutput()); + multiplicationFilter->Update(); + + //FloatImageType::Pointer maskedImage = FloatImageType::New(); + //maskedImage = MultiplicationFilter->GetOutput(); + + FloatWriterType::Pointer floatWriter = FloatWriterType::New(); + std::string s = faFiles.toStdString().append("/"+*outputIt); + floatWriter->SetFileName(s.c_str()); + floatWriter->SetInput(multiplicationFilter->GetOutput()); + floatWriter->Update(); + + ++faIt; + ++maskIt; + ++outputIt; + } +} +/* +void QmitkTractbasedSpatialStatisticsView::Project() +{ + + typedef itk::SkeletonizationFilter SkeletonizationFilter; + SkeletonizationFilter::Pointer skeletonizer = SkeletonizationFilter::New(); + + if (m_CurrentSelection) + { + mitk::DataStorage::SetOfObjects::Pointer set = + mitk::DataStorage::SetOfObjects::New(); + + + mitk::Image::Pointer vol4d; + mitk::TbssImage::Pointer tbssImg; + bool containsSkeleton = false; + bool containsSkeletonMask = false; + bool containsDistanceMap = false; + bool containsGradient = false; + + + for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); + i != m_CurrentSelection->End(); + ++i) + { + if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) + { + mitk::DataNode::Pointer node = nodeObj->GetDataNode(); + + if(QString("Image").compare(node->GetData()->GetNameOfClass())==0) + { + mitk::Image* img = static_cast(node->GetData()); + if(img->GetDimension() == 4) + { + // 4d volume + vol4d = img; + } + } + else if(QString("TbssImage").compare(node->GetData()->GetNameOfClass())==0) + { + mitk::TbssImage::Pointer tempTbss = static_cast(node->GetData()); + if(tempTbss->GetIsMeta()) + { + containsDistanceMap = tempTbss->GetContainsDistanceMap(); + containsSkeleton = tempTbss->GetContainsDistanceMap(); + containsSkeletonMask = tempTbss->GetContainsSkeletonMask(); + containsGradient = tempTbss->GetContainsGradient(); + + + if(containsDistanceMap && containsSkeleton && containsSkeletonMask && containsGradient) + { + tbssImg = tempTbss; + } + } + } + } + } + + if(vol4d.IsNotNull() && tbssImg.IsNotNull()) + { + typedef itk::ProjectionFilter ProjectionFilterType; + ProjectionFilterType::Pointer projector = ProjectionFilterType::New(); + + VectorImageType::Pointer vecImg = ConvertToVectorImage(vol4d); + + + + } + } + + + +} +*/ + + +VectorImageType::Pointer QmitkTractbasedSpatialStatisticsView::ConvertToVectorImage(mitk::Image::Pointer mitkImage) +{ + VectorImageType::Pointer vecImg = VectorImageType::New(); + + mitk::Geometry3D* geo = mitkImage->GetGeometry(); + mitk::Vector3D spacing = geo->GetSpacing(); + mitk::Point3D origin = geo->GetOrigin(); + + VectorImageType::SpacingType vecSpacing; + vecSpacing[0] = spacing[0]; + vecSpacing[1] = spacing[1]; + vecSpacing[2] = spacing[2]; + + VectorImageType::PointType vecOrigin; + vecOrigin[0] = origin[0]; + vecOrigin[1] = origin[1]; + vecOrigin[2] = origin[2]; + + VectorImageType::SizeType size; + size[0] = mitkImage->GetDimension(0); + size[1] = mitkImage->GetDimension(1); + size[2] = mitkImage->GetDimension(2); + + + vecImg->SetSpacing(vecSpacing); + vecImg->SetOrigin(vecOrigin); + vecImg->SetRegions(size); + vecImg->SetVectorLength(mitkImage->GetDimension(3)); + vecImg->Allocate(); + + + for(int x=0; x pixel = vecImg->GetPixel(ix); + + for (int t=0; tGetPixelValueByIndex(ix, t); + pixel.SetElement(t, f); + } + vecImg->SetPixel(ix, pixel); + } + } + } + + return vecImg; + +} +/* +void QmitkTractbasedSpatialStatisticsView::Skeletonize() +{ + + typedef itk::SkeletonizationFilter SkeletonizationFilter; + SkeletonizationFilter::Pointer skeletonizer = SkeletonizationFilter::New(); + + if (m_CurrentSelection) + { + mitk::DataStorage::SetOfObjects::Pointer set = + mitk::DataStorage::SetOfObjects::New(); + + + for (IStructuredSelection::iterator i = m_CurrentSelection->Begin(); + i != m_CurrentSelection->End(); + ++i) + { + + if (mitk::DataNodeObject::Pointer nodeObj = i->Cast()) + { + mitk::DataNode::Pointer node = nodeObj->GetDataNode(); + if(QString("Image").compare(node->GetData()->GetNameOfClass())==0) + { + mitk::Image* img = static_cast(node->GetData()); + if(img->GetDimension() == 3) + { + FloatImageType::Pointer itkImg = FloatImageType::New(); + mitk::CastToItkImage(img, itkImg); + skeletonizer->SetInput(itkImg); + skeletonizer->Update(); + + + // Retrieve skeletonize image and put it in the datamanager + FloatImageType::Pointer skeletonizedImg = skeletonizer->GetOutput(); + mitk::Image::Pointer mitkSkeleton = mitk::Image::New(); + mitk::CastToMitkImage(skeletonizedImg, mitkSkeleton); + AddTbssToDataStorage(mitkSkeleton, "Skeletonized"); + + double threshold = -1.0; + while(threshold == -1.0) + { + // threshold = QInputDialog::getDouble(m_Controls->m_Skeletonize, tr("Specify the FA threshold"), + // tr("Threshold:"), QLineEdit::Normal, + // 0.2); + + if(threshold < 0.0 || threshold > 1.0) + { + QMessageBox msgBox; + msgBox.setText("Please choose a value between 0 and 1"); + msgBox.exec(); + threshold = -1.0; + } + } + + typedef itk::BinaryThresholdImageFilter ThresholdFilterType; + CharImageType::Pointer thresholded = CharImageType::New(); + ThresholdFilterType::Pointer thresholder = ThresholdFilterType::New(); + thresholder->SetInput(skeletonizedImg); + thresholder->SetLowerThreshold(threshold); + thresholder->SetUpperThreshold(std::numeric_limits::max()); + thresholder->SetOutsideValue(0); + thresholder->SetInsideValue(1); + thresholder->Update(); + + + CharImageType::Pointer thresholdedImg = thresholder->GetOutput(); + mitk::Image::Pointer mitkThresholded = mitk::Image::New(); + mitk::CastToMitkImage(thresholdedImg, mitkThresholded); + AddTbssToDataStorage(mitkThresholded, "Skeleton mask"); + + + typedef itk::DistanceMapFilter DistanceMapFilterType; + DistanceMapFilterType::Pointer distanceMapper = DistanceMapFilterType::New(); + distanceMapper->SetInput(thresholdedImg); + distanceMapper->Update(); + + FloatImageType::Pointer distanceMap = distanceMapper->GetOutput(); + mitk::Image::Pointer mitkDistMap = mitk::Image::New(); + mitk::CastToMitkImage(distanceMap, mitkDistMap); + AddTbssToDataStorage(mitkDistMap, "Distance map"); + + DirectionImageType::Pointer dirImage = skeletonizer->GetDirectionImage(); + + mitk::TbssGradientImage::Pointer tbssGradImg = mitk::TbssGradientImage::New(); + tbssGradImg->SetImage(dirImage); + tbssGradImg->InitializeFromVectorImage(); + AddTbssToDataStorage(tbssGradImg, "Gradient image"); + + } + } + } + } + } + + + +} + +/* +void QmitkTractbasedSpatialStatisticsView::InitializeGridByVectorImage() +{ + // Read vector image from file + typedef itk::ImageFileReader< FloatVectorImageType > VectorReaderType; + VectorReaderType::Pointer vectorReader = VectorReaderType::New(); + vectorReader->SetFileName("E:\\tbss\\testing\\Gradient.mhd"); + vectorReader->Update(); + FloatVectorImageType::Pointer directions = vectorReader->GetOutput(); + + + // Read roi from file. + CharReaderType::Pointer roiReader = CharReaderType::New(); + roiReader->SetFileName("E:\\tbss\\testing\\debugging skeletonization\\segment2.mhd"); + roiReader->Update(); + CharImageType::Pointer roi = roiReader->GetOutput(); + + DoInitializeGridByVectorImage(directions, roi, std::string("directions")); + + +} + + + + +void QmitkTractbasedSpatialStatisticsView::DoInitializeGridByVectorImage(FloatVectorImageType::Pointer vectorpic, CharImageType::Pointer roi, std::string name) +{ + //vtkStructuredGrid* grid = vtkStructuredGrid::New(); + itk::Matrix itkdirection = vectorpic->GetDirection(); + itk::Matrix itkinversedirection = itk::Matrix(itkdirection.GetInverse()); + std::vector GridPoints; + vtkPoints *points = vtkPoints::New(); + mitk::Geometry3D::Pointer geom = mitk::Geometry3D::New(); + vtkLinearTransform *vtktransform; + vtkLinearTransform *inverse; + mitk::Image::Pointer geomget = mitk::Image::New(); + geomget->InitializeByItk(vectorpic.GetPointer()); + geom = geomget->GetGeometry(); + vtktransform = geom->GetVtkTransform(); + inverse = vtktransform->GetLinearInverse(); + vtkFloatArray * directions = vtkFloatArray::New(); + directions->SetName("Vectors"); + directions->SetNumberOfComponents(3); + + // Iterator for the vector image + itk::ImageRegionIterator it_input(vectorpic, vectorpic->GetLargestPossibleRegion()); + FloatVectorType nullvector; + nullvector.Fill(0); + double lengthsum = 0; + int id = 0; + + // Iterator for the roi + itk::ImageRegionIterator roiIt(roi, roi->GetLargestPossibleRegion()); + roiIt.GoToBegin(); + + for(it_input.GoToBegin(); !( it_input.IsAtEnd() || roiIt.IsAtEnd() ); ++it_input) + { + //VectorType val = it_input.Value(); + if(it_input.Value() != nullvector && roiIt.Get() != 0) + { + //itk::Point point; + mitk::Point3D mitkpoint, mitkworldpoint; + mitk::Point3D mitkendpoint, mitkworldendpoint; + mitk::Vector3D mitkvector, mitktransvector; + itk::Point direction = it_input.Value().GetDataPointer(); + //itk::Index<3> in_input = it_input.GetIndex(); + //itk::ContinuousIndex cindirection; + FloatVectorType transvec = it_input.Value(); + mitkvector[0] = transvec[0]; + mitkvector[1] = transvec[1]; + mitkvector[2] = transvec[2]; + //mitkvector[2] = 0.0; + + + mitkpoint[0] = it_input.GetIndex()[0]; + mitkpoint[1] = it_input.GetIndex()[1]; + mitkpoint[2] = it_input.GetIndex()[2]; + + mitkendpoint[0] = mitkpoint[0] + mitkvector[0]; + mitkendpoint[1] = mitkpoint[1] + mitkvector[1]; + mitkendpoint[2] = mitkpoint[2] + mitkvector[2]; + + //mitkpoint.setXYZ((ScalarType)point[0],(ScalarType)point[1],(ScalarType)point[2]); + geom->IndexToWorld(mitkpoint, mitkworldpoint); + geom->IndexToWorld(mitkendpoint, mitkworldendpoint); + mitktransvector[0] = mitkworldendpoint[0] - mitkworldpoint[0]; + mitktransvector[1] = mitkworldendpoint[1] - mitkworldpoint[1]; + mitktransvector[2] = mitkworldendpoint[2] - mitkworldpoint[2]; + + lengthsum += mitktransvector.GetNorm(); + + directions->InsertTuple3(id,mitktransvector[0],mitktransvector[1],mitktransvector[2]); + + points->InsertPoint(id,mitkworldpoint[0],mitkworldpoint[1],mitkworldpoint[2]); + id++; + //for (unsigned short loop = 0; (loop < 20) && (!it_input.IsAtEnd()); loop++) + //{ + // ++it_input; + //} + if(it_input.IsAtEnd()) + { + break; + } + } + ++roiIt; + } + double meanlength = lengthsum / id; + vtkGlyph3D* glyph = vtkGlyph3D::New(); + vtkUnstructuredGrid* ugrid = vtkUnstructuredGrid::New(); + ugrid->SetPoints(points); + ugrid->GetPointData()->SetVectors(directions); + glyph->SetInput(ugrid); + glyph->SetScaleModeToScaleByVector(); + glyph->SetScaleFactor(0.5); + glyph->SetColorModeToColorByScalar(); + //glyph->ClampingOn(); + vtkArrowSource* arrow = vtkArrowSource::New(); + if(meanlength > 5) {arrow->SetTipLength(0);arrow->SetTipRadius(0);} + + arrow->SetShaftRadius(0.03/meanlength); + //arrow->SetTipRadius(0.05/meanlength); + + glyph->SetSource(arrow->GetOutput()); + glyph->Update(); + + mitk::Surface::Pointer glyph_surface = mitk::Surface::New(); + + glyph_surface->SetVtkPolyData(glyph->GetOutput()); + glyph_surface->UpdateOutputInformation(); + + mitk::DataNode::Pointer gridNode = mitk::DataNode::New(); + gridNode->SetProperty( "name", mitk::StringProperty::New(name.c_str()) ); + //m_GridNode->SetProperty( "color" , m_GridColor); + gridNode->SetProperty( "visible", mitk::BoolProperty::New(true) ); + gridNode->SetProperty( "segmentation", mitk::BoolProperty::New(true) ); + gridNode->SetProperty( "ID-Tag", mitk::StringProperty::New("grid") ); + gridNode->SetProperty( "shader", mitk::StringProperty::New("mitkShaderLightning") ); + gridNode->SetData( glyph_surface ); + GetDefaultDataStorage()->Add(gridNode); + + +} + +*/ diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h new file mode 100644 index 0000000000..c2eb8883e4 --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h @@ -0,0 +1,252 @@ +/*========================================================================= + +Program: Medical Imaging & Interaction Toolkit +Language: C++ +Date: $Date: 2010-03-31 16:40:27 +0200 (Mi, 31 Mrz 2010) $ +Version: $Revision: 21975 $ + +Copyright (c) German Cancer Research Center, Division of Medical and +Biological Informatics. All rights reserved. +See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ + +#ifndef QmitkTractbasedSpatialStatisticsView_h +#define QmitkTractbasedSpatialStatisticsView_h + + +#include +#include +#include + + +#include + +#include "ui_QmitkTractbasedSpatialStatisticsViewControls.h" +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "QmitkTbssTableModel.h" +#include "QmitkTbssMetaTableModel.h" + + + + +typedef short DiffusionPixelType; +typedef itk::Image CharImageType; +typedef itk::Image UCharImageType; +typedef itk::Image Float4DImageType; +typedef itk::Image FloatImageType; +typedef itk::Vector IntVectorType; +typedef itk::VectorImage DirectionImageType; +typedef itk::VectorImage VectorImageType; + + +typedef itk::ImageFileReader< CharImageType > CharReaderType; +typedef itk::ImageFileReader< UCharImageType > UCharReaderType; +typedef itk::ImageFileWriter< CharImageType > CharWriterType; +typedef itk::ImageFileReader< FloatImageType > FloatReaderType; +typedef itk::ImageFileWriter< FloatImageType > FloatWriterType; +typedef itk::ImageFileReader< Float4DImageType > Float4DReaderType; +typedef itk::ImageFileWriter< Float4DImageType > Float4DWriterType; + + +struct TbssSelListener; + + + +/*! + \brief QmitkTractbasedSpatialStatisticsView + + \warning This application module is not yet documented. Use "svn blame/praise/annotate" and ask the author to provide basic documentation. + + \sa QmitkFunctionalitymitkTbssWorkspaceManager + \ingroup Functionalities +*/ +class QmitkTractbasedSpatialStatisticsView : public QmitkFunctionality +{ + + + friend struct TbssSelListener; + + + + // this is needed for all Qt objesetupUicts that should have a Qt meta-object + // (everything that derives from QObject and wants to have signal/slots) + Q_OBJECT + + public: + + static const std::string VIEW_ID; + + QmitkTractbasedSpatialStatisticsView(); + virtual ~QmitkTractbasedSpatialStatisticsView(); + + virtual void CreateQtPartControl(QWidget *parent); + + /// \brief Creation of the connections of main and control widget + virtual void CreateConnections(); + + virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); + virtual void StdMultiWidgetNotAvailable(); + + /// \brief Called when the functionality is activated + virtual void Activated(); + + virtual void Deactivated(); + + + protected slots: + + + //void OutputValues(); + + + + // void InitializeGridByVectorImage(); + + void Masking(); + + + + void CreateRoi(); + + + void Clustering(); + + void AdjustPlotMeasure(const QString & text); + + void Clicked(const QwtDoublePoint& pos); + + void TbssImport(); + + + void AddGroup(); + void RemoveGroup(); + + void CopyToClipboard(); + + + + protected: + + /// \brief called by QmitkFunctionality when DataManager's selection has changed + virtual void OnSelectionChanged( std::vector nodes ); + + void Plot(mitk::TbssImage*, mitk::TbssRoiImage*); + + void InitPointsets(); + + void SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name); + + berry::ISelectionListener::Pointer m_SelListener; + berry::IStructuredSelection::ConstPointer m_CurrentSelection; + + bool m_IsInitialized; + + mitk::PointSet::Pointer m_PointSetNode; + mitk::DataNode::Pointer m_P1; + + Ui::QmitkTractbasedSpatialStatisticsViewControls* m_Controls; + + QmitkStdMultiWidget* m_MultiWidget; + + + std::vector SortPoints(CharImageType::Pointer roi, CharImageType::IndexType currentPoint); + + bool PointVisited(std::vector points, CharImageType::IndexType point); + + // Modifies the current point by reference and returns true if no more points need to be visited + CharImageType::IndexType FindNextPoint(std::vector pointsVisited, + CharImageType::IndexType currentPoint, CharImageType::Pointer roi, bool &ready); + + + + + //void DoInitializeGridByVectorImage(FloatVectorImageType::Pointer vectorpic, CharImageType::Pointer roi ,std::string name); + + + + // Tokenizer needed for the roi files + void Tokenize(const std::string& str, + std::vector& tokens, + const std::string& delimiters = " ") + { + // Skip delimiters at beginning. + std::string::size_type lastPos = str.find_first_not_of(delimiters, 0); + // Find first "non-delimiter". + std::string::size_type pos = str.find_first_of(delimiters, lastPos); + + while (std::string::npos != pos || std::string::npos != lastPos) + { + // Found a token, add it to the vector. + tokens.push_back(str.substr(lastPos, pos - lastPos)); + // Skip delimiters. Note the "not_of" + lastPos = str.find_first_not_of(delimiters, pos); + // Find next "non-delimiter" + pos = str.find_first_of(delimiters, lastPos); + } + } + + mitk::DataNode::Pointer readNode(std::string f) + { + mitk::DataNode::Pointer node; + mitk::DataNodeFactory::Pointer nodeReader = mitk::DataNodeFactory::New(); + try + { + nodeReader->SetFileName(f); + nodeReader->Update(); + node = nodeReader->GetOutput(); + } + catch(...) { + MITK_ERROR << "Could not read file"; + return NULL; + } + + return node; + } + + /*template < typename TPixel, unsigned int VImageDimension > + void ToITK4D( itk::Image* inputImage, Float4DImageType::Pointer& outputImage );*/ + + + std::string ReadFile(std::string whatfile); + + + std::vector< itk::Index<3> > m_Roi; + + std::string m_CurrentStructure; + + mitk::Geometry3D* m_CurrentGeometry; + + QmitkTbssTableModel* m_GroupModel; + + + + void AddTbssToDataStorage(mitk::Image* image, std::string name); + + mitk::TbssImage::Pointer m_CurrentTbssMetaImage; + + VectorImageType::Pointer ConvertToVectorImage(mitk::Image::Pointer mitkImg); + +}; + + + +#endif // _QMITKTRACTBASEDSPATIALSTATISTICSVIEW_H_INCLUDED + diff --git a/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsViewControls.ui b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsViewControls.ui new file mode 100644 index 0000000000..0ac2df05ff --- /dev/null +++ b/Modules/Bundles/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsViewControls.ui @@ -0,0 +1,402 @@ + + + QmitkTractbasedSpatialStatisticsViewControls + + + + 0 + 0 + 431 + 811 + + + + + 0 + 0 + + + + QmitkTemplate + + + + + + FSL import + + + false + + + false + + + + + + Here subject data and tbss meta data can be imported from FSL into the MITK TBSS module + + + 0 + + + + Subject data + + + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + + + Group information + + + QAbstractItemView::SelectRows + + + + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + + + Add group entry. After that give the group a name and the correct number + + + Add + + + + + + + Remove selected entries + + + Remove + + + + + + + false + + + Import a 4D image containing group data after group information has been set + + + Import subject data + + + + + + + + + + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + + + Diffusion measure + + + + + + + Measurement in the to be imported 4D image + + + Fractional Anisotropy + + + + + + + + + + + Meta data + + + + + + true + + + Meta data + + + QAbstractItemView::SelectRows + + + + + + + QFrame::StyledPanel + + + QFrame::Raised + + + + + + Select the function of the tbss image + + + + Add skeleton mask + + + + + Add mean FA Skeleton + + + + + Add gradient image + + + + + Add tubular structure + + + + + Add distance map + + + + + + + + Remove meta data entries + + + Remove + + + + + + + Create a tbss meta image with the images in the table. + + + Create TBSS meta file + + + + + + + + + + + + + + + + + Tract-specific analysis + + + + + + true + + + To create a roi first load a tbss meta image into the datamanager + + + 1 + + + + ROIs + + + + QFormLayout::AllNonFixedFieldsGrow + + + + + current selection + + + mean FA skeleton: + + + + + + + Points on Roi + + + + + + + + 100 + 100 + + + + + 0 + 100 + + + + Use this widget to create points on the ROI by shift-leftclick on the right positions on the skeleton. Then click Create Roi. The Roi that will be created will pass through the points in the order of occurence in this list + + + + + + + false + + + No suitable tbss meta image selected yet. The meta image needs to contain a mean FA skeleton and a skeleton mask + + + Create ROI + + + + + + + + 0 + 0 + + + + Points on the ROI + + + + + + + Name + + + + + + + Give a name to the region of interest + + + roiname + + + + + + + Structure info + + + + + + + On what anatomical structure lies the ROI? + + + Structure + + + + + + + + Measuring + + + + + + To plot, load a tbss image with subject information and a region of interest corresponding to the study and select them both + + + + + + + Copy to clipboard + + + + + + + + + + + + + + + + QmitkPointListWidget + QWidget +
QmitkPointListWidget.h
+ + + QmitkTbssRoiAnalysisWidget + QWidget +
QmitkTbssRoiAnalysisWidget.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 5ff89f0889..988aa100a2 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,52 +1,54 @@ #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/QmitkGibbsTrackingView.h" #include "src/internal/QmitkStochasticFiberTrackingView.h" #include "src/internal/QmitkFiberProcessingView.h" #include "src/internal/QmitkFiberBundleDeveloperView.h" #include "src/internal/QmitkPartialVolumeAnalysisView.h" #include "src/internal/QmitkIVIMView.h" #include "src/internal/QmitkScreenshotMaker.h" +#include "src/internal/QmitkTractbasedSpatialStatisticsView.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(QmitkGibbsTrackingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkStochasticFiberTrackingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkFiberProcessingView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkFiberBundleDeveloperView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkPartialVolumeAnalysisView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkIVIMView, context) BERRY_REGISTER_EXTENSION_CLASS(QmitkScreenshotMaker, context) + BERRY_REGISTER_EXTENSION_CLASS(QmitkTractbasedSpatialStatisticsView, context) } void PluginActivator::stop(ctkPluginContext* context) { Q_UNUSED(context) } } Q_EXPORT_PLUGIN2(org_mitk_gui_qt_diffusionimaging, mitk::PluginActivator) diff --git a/Modules/CMakeLists.txt b/Modules/CMakeLists.txt index 569ca82327..e062f6c3d6 100644 --- a/Modules/CMakeLists.txt +++ b/Modules/CMakeLists.txt @@ -1,46 +1,47 @@ SET(LIBPOSTFIX "Ext") SET(module_dirs SceneSerializationBase PlanarFigure ImageExtraction ImageStatistics MitkExt SceneSerialization QmitkExt + GraphAlgorithms DiffusionImaging GPGPU IGT CameraCalibration IGTUI RigidRegistration RigidRegistrationUI DeformableRegistration DeformableRegistrationUI OpenCVVideoSupport Overlays InputDevices ToFHardware ToFProcessing - ToFUI + ToFUI ) SET(MITK_DEFAULT_SUBPROJECTS MITK-Modules) FOREACH(module_dir ${module_dirs}) ADD_SUBDIRECTORY(${module_dir}) ENDFOREACH() IF(MITK_PRIVATE_MODULES) FILE(GLOB all_subdirs RELATIVE ${MITK_PRIVATE_MODULES} ${MITK_PRIVATE_MODULES}/*) FOREACH(subdir ${all_subdirs}) STRING(FIND ${subdir} "." _result) IF(_result EQUAL -1) IF(EXISTS ${MITK_PRIVATE_MODULES}/${subdir}/CMakeLists.txt) MESSAGE(STATUS "Found private module ${subdir}") ADD_SUBDIRECTORY(${MITK_PRIVATE_MODULES}/${subdir} private_modules/${subdir}) ENDIF() ENDIF() ENDFOREACH() ENDIF(MITK_PRIVATE_MODULES) diff --git a/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.cpp b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.cpp new file mode 100644 index 0000000000..d5e9be0f60 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.cpp @@ -0,0 +1,162 @@ + +#ifndef __mitkTractAnalyzer_cpp +#define __mitkTractAnalyzer_cpp + + + +#include +#include + +#include + + +#include +#include + +#include +#include + + +using namespace std; + +namespace mitk { + + TractAnalyzer::TractAnalyzer() + { + + } + + + void TractAnalyzer::BuildGraph(itk::Index<3> startPoint, itk::Index<3> endPoint) + { + + typedef itk::ShortestPathImageFilter ShortestPathFilterType; + typedef itk::ShortestPathCostFunctionTbss CostFunctionType; + + + FloatImageType::Pointer meanSkeleton; + + mitk::CastToItkImage(m_InputImage, meanSkeleton); + + // Only use the mitk image + + + + if(meanSkeleton) + { + CostFunctionType::Pointer costFunction = CostFunctionType::New(); + costFunction->SetImage(meanSkeleton); + costFunction->SetStartIndex(startPoint); + costFunction->SetEndIndex(endPoint); + costFunction->SetThreshold(m_Threshold); + + ShortestPathFilterType::Pointer pathFinder = ShortestPathFilterType::New(); + pathFinder->SetCostFunction(costFunction); + pathFinder->SetFullNeighborsMode(true); + //pathFinder->SetCalcMode(ShortestPathFilterType::A_STAR); + pathFinder->SetInput(meanSkeleton); + pathFinder->SetStartIndex(startPoint); + pathFinder->SetEndIndex(endPoint); + pathFinder->Update(); + + m_Path = pathFinder->GetVectorPath(); + + + m_RoiImg = pathFinder->GetOutput(); + } + + + + } + + void TractAnalyzer::MeasureRoi() + { + + // Output two types + ProjectionsImageType::SizeType size = m_Projections->GetLargestPossibleRegion().GetSize(); + + std::ofstream file(m_FileName.c_str()); + + std::vector individuals; + for(int i=0; i group = m_Groups[i]; + for(int j=0; j ix = m_Roi[k]; + itk::Index<4> ix4; + ix4[0] = ix[0]; + ix4[1] = ix[1]; + ix4[2] = ix[2]; + ix4[3] = j; + + float f = m_Projections->GetPixel(ix4); + + file << f << " "; + + } + + file << "\n"; + + } + + file.close(); + + + + // Write the long format output + std::ofstream fileLong(m_FileNameLong.c_str()); + + fileLong << "ID " << "group " << "position " << "value\n"; + + for(int i=0; i ix = m_Roi[i]; + itk::Index<4> ix4; + ix4[0] = ix[0]; + ix4[1] = ix[1]; + ix4[2] = ix[2]; + ix4[3] = j; + + float f = m_Projections->GetPixel(ix4); + + fileLong << "ID" << j << " " << individuals[j] << " " << "pos" << i << " "<< f << "\n"; + } + } + + fileLong.close(); + + + } + + + +} +#endif diff --git a/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.h b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.h new file mode 100644 index 0000000000..6fbc00d524 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/mitkTractAnalyzer.h @@ -0,0 +1,144 @@ +/*========================================================================= + + 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 __mitkTractAnalyzer_h_ +#define __mitkTractAnalyzer_h_ + +#include "MitkDiffusionImagingMBIExports.h" +#include +#include + + +namespace mitk{ +/** \class TractAnalyzer + */ + + +class MitkDiffusionImagingMBI_EXPORT TractAnalyzer +{ + +public: + + + TractAnalyzer(); + ~TractAnalyzer() {}; + + typedef itk::Image CharImageType; + typedef itk::Image FloatImageType; + typedef itk::Image ProjectionsImageType; + typedef itk::VectorImage VectorImageType; + + + + /* + void SetSkeleton(CharImageType::Pointer skeleton) + { + m_Skeleton = skeleton; + } + + void SetMeanSkeleton(FloatImageType::Pointer i) + { + m_MeanSkeleton = i; + }*/ + + + void SetTbssImage(mitk::TbssImage::Pointer tbssImg) + { + m_TbssImage = tbssImg; + } + + void SetProjections(ProjectionsImageType::Pointer projections) + { + m_Projections = projections; + } + + void BuildGraph(itk::Index<3> startPoint, itk::Index<3> endPoint); + + std::vector< itk::Index<3> > GetPath() + { + return m_Path; + } + + void SetFileName(std::string fname) + { + m_FileName = fname; + } + + void SetFileNameLong(std::string fname) + { + m_FileNameLong = fname; + } + + void SetRoi(std::vector< itk::Index<3> > roi) + { + m_Roi = roi; + } + + CharImageType::Pointer GetRoiImage() + { + return m_RoiImg; + } + + void SetGroups(std::vector< std::pair > groups) + { + m_Groups = groups; + } + + void MeasureRoi(); + + void SetInputImage(mitk::Image::Pointer inputImage) + { + m_InputImage = inputImage; + } + + + void SetThreshold(double threshold) + { + m_Threshold = threshold; + } + + +protected: + + + //CharImageType::Pointer m_Skeleton; + CharImageType::Pointer m_RoiImg; + ProjectionsImageType::Pointer m_Projections; + //FloatImageType::Pointer m_MeanSkeleton; + mitk::TbssImage::Pointer m_TbssImage; + + mitk::Image::Pointer m_InputImage; + + double m_Threshold; + + std::vector< itk::Index<3> > m_Path; + + std::string m_FileName; + + std::string m_FileNameLong; // For the regression analysis 'friendly' file + + std::vector< std::pair > m_Groups; + + std::vector< itk::Index<3> > m_Roi; + +private: + +}; + +} + +#endif //__itkTractAnalyzer_h_ + diff --git a/Modules/DiffusionImaging/CMakeLists.txt b/Modules/DiffusionImaging/CMakeLists.txt index 5c587e23ce..c5e4160b60 100644 --- a/Modules/DiffusionImaging/CMakeLists.txt +++ b/Modules/DiffusionImaging/CMakeLists.txt @@ -1,28 +1,28 @@ FIND_PACKAGE(ITK) IF(ITK_GDCM_DIR) INCLUDE(${ITK_GDCM_DIR}/GDCMConfig.cmake) IF(GDCM_MAJOR_VERSION EQUAL 2) ADD_DEFINITIONS(-DGDCM2) SET(ITK_USES_GDCM2 1) ENDIF(GDCM_MAJOR_VERSION EQUAL 2) ENDIF(ITK_GDCM_DIR) MITK_CREATE_MODULE( MitkDiffusionImaging SUBPROJECTS MITK-DTI INCLUDE_DIRS IODataStructures Reconstruction Tractography Rendering Algorithms DicomImport Interactions IODataStructures/DiffusionWeightedImages IODataStructures/QBallImages IODataStructures/TensorImages IODataStructures/FiberBundle IODataStructures/FiberBundleX IODataStructures/PlanarFigureComposite IODataStructures/TbssImages ${CMAKE_CURRENT_BINARY_DIR} - DEPENDS MitkExt SceneSerializationBase QmitkExt + DEPENDS MitkExt SceneSerializationBase QmitkExt MitkGraphAlgorithms PACKAGE_DEPENDS Boost ) MITK_USE_MODULE(MitkDiffusionImaging) if(MitkDiffusionImaging_IS_ENABLED) file(DOWNLOAD http://mitk.org/download/data/FibertrackingLUT.tar.gz ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/FibertrackingLUT.tar.gz TIMEOUT 10) execute_process(COMMAND ${CMAKE_COMMAND} -E chdir ${CMAKE_RUNTIME_OUTPUT_DIRECTORY} tar xzf FibertrackingLUT.tar.gz) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/Rendering/mitkShaderFiberClipping.xml ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/mitkShaderFiberClipping.xml) MITK_INSTALL(FILES ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/mitkShaderFiberClipping.xml ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/FiberTrackingLUTBaryCoords.bin ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/FiberTrackingLUTIndices.bin) endif() ADD_SUBDIRECTORY(Testing) CONFIGURE_FILE(${CMAKE_CURRENT_SOURCE_DIR}/mitkDiffusionImagingConfigure.h.in ${CMAKE_CURRENT_BINARY_DIR}/mitkDiffusionImagingConfigure.h) diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 961349c228..8f7c4358de 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,180 +1,181 @@ SET(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundle IODataStructures/FiberBundle/mitkFiberBundle.cpp IODataStructures/FiberBundle/mitkFiberBundleWriter.cpp IODataStructures/FiberBundle/mitkFiberBundleReader.cpp IODataStructures/FiberBundle/mitkFiberBundleIOFactory.cpp IODataStructures/FiberBundle/mitkFiberBundleWriterFactory.cpp IODataStructures/FiberBundle/mitkFiberBundleSerializer.cpp IODataStructures/FiberBundle/mitkParticle.cpp IODataStructures/FiberBundle/mitkParticleGrid.cpp # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures -> Tbss IODataStructures/TbssImages/mitkTbssImageSource.cpp IODataStructures/TbssImages/mitkTbssRoiImageSource.cpp IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageIOFactory.cpp IODataStructures/TbssImages/mitkTbssImage.cpp IODataStructures/TbssImages/mitkTbssRoiImage.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriterFactory.cpp IODataStructures/TbssImages/mitkTbssImporter.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleMapper2D.cpp Rendering/mitkFiberBundleMapper3D.cpp Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp Rendering/mitkTbssImageMapper.cpp Rendering/mitkPlanarCircleMapper3D.cpp Rendering/mitkPlanarPolygonMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp + Algorithms/mitkTractAnalyzer.cpp # Tractography Tractography/itkStochasticTractographyFilter.h ) SET(H_FILES # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkTbssImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleMapper2D.h Rendering/mitkFiberBundleMapper3D.h Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h Rendering/mitkPlanarCircleMapper3D.h Rendering/mitkPlanarPolygonMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h Reconstruction/itkRegularizedIVIMReconstructionFilter.h Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/FiberBundle/itkSlowPolyLineParametricPath.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h # Tractography Tractography/itkGibbsTrackingFilter.h Tractography/itkStochasticTractographyFilter.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/itkTractsToFiberEndingsImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkSkeletonizationFilter.h ) SET( TOOL_FILES ) IF(WIN32) ENDIF(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp ) diff --git a/Modules/GraphAlgorithms/CMakeLists.txt b/Modules/GraphAlgorithms/CMakeLists.txt new file mode 100644 index 0000000000..efdc08f77a --- /dev/null +++ b/Modules/GraphAlgorithms/CMakeLists.txt @@ -0,0 +1,3 @@ +# CREATE MODULE HERE +MITK_CREATE_MODULE(MitkGraphAlgorithms +DEPENDS Mitk ImageStatistics ) diff --git a/Modules/GraphAlgorithms/files.cmake b/Modules/GraphAlgorithms/files.cmake new file mode 100644 index 0000000000..c4c60e7126 --- /dev/null +++ b/Modules/GraphAlgorithms/files.cmake @@ -0,0 +1,11 @@ +SET(CPP_FILES + itkShortestPathNode.cpp +) + +SET(H_FILES + itkShortestPathCostFunction.h + itkShortestPathCostFunctionTbss.h + itkShortestPathNode.h + itkShortestPathImageFilter.h + +) diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunction.h b/Modules/GraphAlgorithms/itkShortestPathCostFunction.h new file mode 100644 index 0000000000..a16d260bee --- /dev/null +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunction.h @@ -0,0 +1,85 @@ +#ifndef __itkShortestPathCostFunction_h +#define __itkShortestPathCostFunction_h + +#include "itkObject.h" +#include "itkObjectFactory.h" +#include "itkShapedNeighborhoodIterator.h" +#include + +#include + + +namespace itk +{ + + // \brief this is a pure virtual superclass for all cost function for the itkShortestPathImageFilter + template + class ITK_EXPORT ShortestPathCostFunction : public Object + { + public: + + /** Standard class typedefs. */ + typedef ShortestPathCostFunction Self; + typedef Object Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef ShapedNeighborhoodIterator< TInputImageType > itkShapedNeighborhoodIteratorType; + + /** Run-time type information (and related methods). */ + itkTypeMacro(ShortestPathCostFunction, Object); + + /** Type definition for the input image. */ + typedef TInputImageType ImageType; + + // More typdefs for convenience + typedef typename TInputImageType::Pointer ImagePointer; + typedef typename TInputImageType::ConstPointer ImageConstPointer; + + + typedef typename TInputImageType::PixelType PixelType; + + typedef typename TInputImageType::IndexType IndexType; + + /** Set the input image. */ + itkSetConstObjectMacro(Image,TInputImageType); + + // \brief calculates the costs for going from pixel1 to pixel2 + virtual double GetCost( IndexType p1, IndexType p2) = 0; + + // \brief returns the minimal costs possible (needed for A*) + virtual double GetMinCost() = 0; + + // \brief Initialize the metric + virtual void Initialize () = 0; + + // \brief Set Starpoint for Path + void SetStartIndex (const IndexType & StartIndex); + + // \brief Set Endpoint for Path + void SetEndIndex(const IndexType & EndIndex); + + + ShortestPathCostFunction(); + + protected: + + virtual ~ShortestPathCostFunction() {}; + void PrintSelf(std::ostream& os, Indent indent) const; + ImageConstPointer m_Image; + IndexType m_StartIndex, m_EndIndex; + + private: + + ShortestPathCostFunction(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + }; + +} // end namespace itk + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkShortestPathCostFunction.txx" +#endif + +#endif /* __itkShortestPathCostFunction_h */ diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunction.txx b/Modules/GraphAlgorithms/itkShortestPathCostFunction.txx new file mode 100644 index 0000000000..f3e9896bd9 --- /dev/null +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunction.txx @@ -0,0 +1,51 @@ +#ifndef __itkShortestPathCostFunction_txx +#define __itkShortestPathCostFunction_txx + +#include "itkShortestPathCostFunction.h" + + +namespace itk +{ + + // Constructor + template + ShortestPathCostFunction + ::ShortestPathCostFunction() + { + } + + + template + void + ShortestPathCostFunction + ::PrintSelf( std::ostream& os, Indent indent ) const + { + Superclass::PrintSelf(os,indent); + } + + + template + void ShortestPathCostFunction:: + SetStartIndex (const typename TInputImageType::IndexType &StartIndex) + { + for (unsigned int i=0;i + void ShortestPathCostFunction:: + SetEndIndex (const typename TInputImageType::IndexType &EndIndex) + { + for (unsigned int i=0;i + + + +namespace itk +{ + + template + class ITK_EXPORT ShortestPathCostFunctionTbss : public ShortestPathCostFunction + { + public: + /** Standard class typedefs. */ + typedef ShortestPathCostFunctionTbss Self; + typedef ShortestPathCostFunction Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + typedef itk::ImageRegionConstIterator ConstIteratorType; + typedef typename TInputImageType::IndexType IndexType; + + + typedef itk::Image FloatImageType; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro(ShortestPathCostFunctionTbss, Object); + + // \brief calculates the costs for going from p1 to p2 + virtual double GetCost( IndexType p1, IndexType p2); + + // \brief Initialize the metric + virtual void Initialize (); + + // \brief returns the minimal costs possible (needed for A*) + virtual double GetMinCost(); + + + ShortestPathCostFunctionTbss(); + + void SetThreshold(double t) + { + m_Threshold = t; + } + + + protected: + + virtual ~ShortestPathCostFunctionTbss() {}; + + double m_Threshold; + + + private: + + + + + }; + +} // end namespace itk + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkShortestPathCostFunctionTbss.txx" +#endif + +#endif /* __itkShortestPathCostFunctionTbss_h */ diff --git a/Modules/GraphAlgorithms/itkShortestPathCostFunctionTbss.txx b/Modules/GraphAlgorithms/itkShortestPathCostFunctionTbss.txx new file mode 100644 index 0000000000..41e6df7fab --- /dev/null +++ b/Modules/GraphAlgorithms/itkShortestPathCostFunctionTbss.txx @@ -0,0 +1,74 @@ +#ifndef __itkShortestPathCostFunctionTbss_txx +#define __itkShortestPathCostFunctionTbss_txx + +#include "itkShortestPathCostFunctionTbss.h" +#include + + + +namespace itk +{ + + // Constructor + template + ShortestPathCostFunctionTbss + ::ShortestPathCostFunctionTbss() + { + } + + + template + double ShortestPathCostFunctionTbss + ::GetCost(IndexType p1 ,IndexType p2) + { + // Very simple Metric: + // The squared difference of both values is defined as cost + + double a,b,c; + ConstIteratorType it( this->m_Image, this->m_Image->GetRequestedRegion()); + + it.SetIndex(p1); + a = it.Get(); + it.SetIndex(p2); + b = it.Get(); + + + if(this->m_Image->GetPixel(p2) < m_Threshold) + { + c = std::numeric_limits::max( ); + } + else + { + double dxSqt = (p1[0]-p2[0]) * (p1[0]-p2[0]);// * (p1[0]-p2[0]); + double dySqt = (p1[1]-p2[1]) * (p1[1]-p2[1]); + double dzSqt = (p1[2]-p2[2]) * (p1[2]-p2[2]); + + double weight = this->m_Image->GetPixel(p2); + + c = (dxSqt + dySqt + dzSqt) + 10000 * (1-weight); + } + + + return c; + } + + + template + void ShortestPathCostFunctionTbss + ::Initialize() + { + } + + template + double ShortestPathCostFunctionTbss + ::GetMinCost() + { + return 1; + } + + + + +} // end namespace itk + +#endif // __itkShortestPathCostFunctionSimple_txx diff --git a/Modules/GraphAlgorithms/itkShortestPathImageFilter.h b/Modules/GraphAlgorithms/itkShortestPathImageFilter.h new file mode 100644 index 0000000000..e4f839813a --- /dev/null +++ b/Modules/GraphAlgorithms/itkShortestPathImageFilter.h @@ -0,0 +1,215 @@ +#ifndef __itkShortestPathImageFilter_h_ +#define __itkShortestPathImageFilter_h_ + +#include "itkImageToImageFilter.h" +#include "itkShortestPathCostFunction.h" +#include "itkShortestPathNode.h" +#include + +#include + +// ------- INFORMATION ---------- +/// SET FUNCTIONS +//void SetInput( ItkImage ) // Compulsory +//void SetStartIndex (const IndexType & StartIndex); // Compulsory +//void SetEndIndex(const IndexType & EndIndex); // Compulsory +//void SetFullNeighborsMode(bool) // Optional (default=false), if false N4, if true N26 +//void SetActivateTimeOut(bool) // Optional (default=false), for debug issues: after 30s algorithms terminates. You can have a look at the VectorOrderImage to see how far it came +//void SetMakeOutputImage(bool) // Optional (default=true), Generate an outputimage of the path. You can also get the path directoy with GetVectorPath() +//void SetCalcAllDistances(bool) // Optional (default=false), Calculate Distances over the whole image. CAREFUL, algorithm time extends a lot. Necessary for GetDistanceImage +//void SetStoreVectorOrder(bool) // Optional (default=false), Stores in which order the pixels were checked. Necessary for GetVectorOrderImage +//void AddEndIndex(const IndexType & EndIndex) //Optional. By calling this function you can add several endpoints! The algorithm will look for several shortest Pathes. From Start to all Endpoints. +// +/// GET FUNCTIONS +//std::vector< itk::Index<3> > GetVectorPath(); // returns the shortest path as vector +//std::vector< std::vector< itk::Index<3> > GetMultipleVectorPathe(); // returns a vector of shortest Pathes (which are vectors of points) +//GetDistanceImage // Returns the distance image +//GetVectorOrderIMage // Returns the Vector Order image +// +// EXAMPLE USE +// pleae see qmitkmitralvalvesegmentation4dtee bundle + + + + + +namespace itk +{ + + template + class ITK_EXPORT ShortestPathImageFilter : + public ImageToImageFilter + { + public: + //Standard Typedefs + typedef ShortestPathImageFilter Self; + typedef ImageToImageFilter Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + // Typdefs for metric + typedef ShortestPathCostFunction< TInputImageType > CostFunctionType; + typedef typename CostFunctionType::Pointer CostFunctionTypePointer; + + + // More typdefs for convenience + typedef TInputImageType InputImageType; + typedef typename TInputImageType::Pointer InputImagePointer; + typedef typename TInputImageType::PixelType InputImagePixelType; + typedef typename TInputImageType::SizeType InputImageSizeType; + typedef typename TInputImageType::IndexType IndexType; + typedef typename itk::ImageRegionIteratorWithIndex< InputImageType > InputImageIteratorType; + + typedef TOutputImageType OutputImageType; + typedef typename TOutputImageType::Pointer OutputImagePointer; + typedef typename TOutputImageType::PixelType OutputImagePixelType; + typedef typename TOutputImageType::IndexType OutputImageIndexType; + typedef ImageRegionIteratorWithIndex< OutputImageType > OutputImageIteratorType; + typedef itk::ShapedNeighborhoodIterator< TInputImageType > itkShapedNeighborhoodIteratorType; + + + // New Macro for smartpointer instantiation + itkNewMacro(Self); + + // Run-time type information + itkTypeMacro(ShortestPathImageFilter, ImageToImageFilter); + + // Display + void PrintSelf( std::ostream& os, Indent indent ) const; + + // Compare function for A_STAR + struct CompareNodeStar + { + bool operator()(ShortestPathNode *a, ShortestPathNode *b) + { + return (a->distAndEst > b->distAndEst); + } + }; + + // \brief Set Starpoint for ShortestPath Calculation + void SetStartIndex (const IndexType & StartIndex); + + // \brief Adds Endpoint for multiple ShortestPath Calculation + void AddEndIndex (const IndexType & index); + + // \brief Set Endpoint for ShortestPath Calculation + void SetEndIndex(const IndexType & EndIndex); + + + // \brief Set FullNeighborsMode. false = no diagonal neighbors, in 2D this means N4 Neigborhood. true = would be N8 in 2D + itkSetMacro (FullNeighborsMode, bool); + itkGetMacro (FullNeighborsMode, bool); + + // \brief (default=true), Produce output image, which shows the shortest path. But you can also get the shortest Path directly as vector with the function GetVectorPath + itkSetMacro (MakeOutputImage, bool); + itkGetMacro (MakeOutputImage, bool); + + // \brief (default=false), Store an Vector of Order, so you can call getVectorOrderImage after update + itkSetMacro (StoreVectorOrder, bool); + itkGetMacro (StoreVectorOrder, bool); + + // \brief (default=false), // Calculate all Distances to all pixels, so you can call getDistanceImage after update (warning algo will take a long time) + itkSetMacro (CalcAllDistances, bool); + itkGetMacro (CalcAllDistances, bool); + + // \brief (default=false), for debug issues: after 30s algorithms terminates. You can have a look at the VectorOrderImage to see how far it came + itkSetMacro (ActivateTimeOut, bool); + itkGetMacro (ActivateTimeOut, bool); + + // \brief returns shortest Path as vector + std::vector< itk::Index<3> > GetVectorPath(); + + // \brief returns Multiple shortest Paths. You can call this function, when u performed a multiple shortest path search (one start, several ends) + std::vector< std::vector< itk::Index<3> > > GetMultipleVectorPaths(); + + // \brief returns the vector order image. It shows in which order the pixels were checked. good for debugging. Be sure to have m_StoreVectorOrder=true + OutputImagePointer GetVectorOrderImage(); + + // \brief returns the distance image. It shows the distances from the startpoint to all other pixels. Be sure to have m_CalcAllDistances=true + OutputImagePointer GetDistanceImage(); + + // \brief cleans up the filter + void CleanUp(); + + itkSetObjectMacro( CostFunction, CostFunctionType ); // itkSetObjectMacro = set function that uses pointer as parameter + itkGetObjectMacro( CostFunction, CostFunctionType ); + + protected: + + std::vector< IndexType > m_endPoints; // if you fill this vector, the algo will not rest until all endPoints have been reached + std::vector< IndexType > m_endPointsClosed; + + ShortestPathNode* m_Nodes; // main list that contains all nodes + NodeNumType m_Graph_NumberOfNodes; + NodeNumType m_Graph_StartNode; + NodeNumType m_Graph_EndNode; + int m_ImageDimensions; + bool m_Graph_fullNeighbors; + std::vector m_Graph_DiscoveredNodeList; + ShortestPathImageFilter(Self&); // intentionally not implemented + void operator=(const Self&); // intentionally not implemented + const static int BACKGROUND = 0; + const static int FOREGROUND = 255; + bool m_FullNeighborsMode; + + bool m_MakeOutputImage; + bool m_StoreVectorOrder; // Store an Vector of Order, so you can call getVectorOrderImage after update + bool m_CalcAllDistances; // Calculate all Distances, so you can call getDistanceImage after update (warning algo will take a long time) + + bool multipleEndPoints; + + bool m_ActivateTimeOut; // if true, then i search max. 30 secs. then abort + CostFunctionTypePointer m_CostFunction; + IndexType m_StartIndex, m_EndIndex; + std::vector< itk::Index<3> > m_VectorPath; + std::vector< std::vector< itk::Index<3> > > m_MultipleVectorPaths; + + + std::vector< NodeNumType > m_VectorOrder; + + + + ShortestPathImageFilter (); + + // \brief Fill m_VectorPath + void MakeShortestPathVector(); + + // \brief Create all the outputs + void MakeOutputs(); + + // \brief Generate Data + void GenerateData(); + + // \brief gets the estimate costs from pixel a to target. + double getEstimatedCostsToTarget(const IndexType & a); + + typename InputImageType::Pointer m_magnitudeImage; + + // \brief Convert a indexnumber of a node in m_Nodes to image coordinates + typename TInputImageType::IndexType NodeToCoord(NodeNumType); + + // \brief Convert image coordinate to a indexnumber of a node in m_Nodes + unsigned int CoordToNode(IndexType); + + // \brief Returns the neighbors of a node + std::vector GetNeighbors(NodeNumType nodeNum, bool FullNeighbors); + + // \brief Check if coords are in bounds of image + bool CoordIsInBounds(IndexType); + + // \brief Initializes the graph + void InitGraph(); + + // \brief Start ShortestPathSearch + void StartShortestPathSearch(); + + }; + +} // end of namespace itk + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkShortestPathImageFilter.txx" +#endif + +#endif diff --git a/Modules/GraphAlgorithms/itkShortestPathImageFilter.txx b/Modules/GraphAlgorithms/itkShortestPathImageFilter.txx new file mode 100644 index 0000000000..14dc080536 --- /dev/null +++ b/Modules/GraphAlgorithms/itkShortestPathImageFilter.txx @@ -0,0 +1,917 @@ +#include "itkShortestPathImageFilter.h" +#include "time.h" +#include "mitkMemoryUtilities.h" +#include +#include +#include + +namespace itk +{ + // Constructor (initialize standard values) + template + ShortestPathImageFilter + ::ShortestPathImageFilter() : + m_FullNeighborsMode(false), + m_MakeOutputImage(true), + m_StoreVectorOrder(false), + m_CalcAllDistances(false), + m_ActivateTimeOut(false), + multipleEndPoints(false), + m_Nodes(0), + m_Graph_NumberOfNodes(0) + { + m_endPoints.clear(); + m_endPointsClosed.clear(); + + if (m_MakeOutputImage) + { + this->SetNumberOfRequiredOutputs(1); + this->SetNthOutput( 0, OutputImageType::New() ); + } + } + + template + inline typename ShortestPathImageFilter::IndexType + ShortestPathImageFilter + ::NodeToCoord (NodeNumType node) + { + const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize(); + int dim = InputImageType::ImageDimension; + IndexType coord; + if (dim == 2) + { + coord[1] = node / size[0]; + coord[0] = node % size[0]; + if ((coord[0] >= size[0]) || (coord[1] >= size[1])) + { + coord[0] = 0; + coord[1] = 0; + } + } + if (dim == 3) + { + coord[2] = node / (size[0]*size[1]); + coord[1] = (node % (size[0]*size[1])) / size[0]; + coord[0] = (node % (size[0]*size[1])) % size[0]; + if ((coord[0] >= size[0]) || (coord[1] >= size[1]) || (coord[2] >= size[2])) + { + coord[0] = 0; + coord[1] = 0; + coord[2] = 0; + } + } + + return coord; + } + + template + inline typename itk::NodeNumType + ShortestPathImageFilter:: + CoordToNode (IndexType coord) + { + const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize(); + int dim = InputImageType::ImageDimension; + NodeNumType node = 0; + if (dim == 2) + { + node = (coord[1]*size[0]) + coord[0]; + } + if (dim == 3) + { + node = (coord[2]*size[0]*size[1]) + (coord[1]*size[0]) + coord[0]; + } + if ((m_Graph_NumberOfNodes > 0) && (node >= m_Graph_NumberOfNodes)) + { + MITK_INFO << "WARNING! Coordinates outside image!"; + MITK_INFO << "Coords = " << coord ; + MITK_INFO << "ImageDim = " << dim ; + MITK_INFO << "RequestedRegionSize = " << size ; + node = 0; + } + + return node; + } + + template + inline bool + ShortestPathImageFilter:: + CoordIsInBounds (IndexType coord) + { + const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize(); + int dim = InputImageType::ImageDimension; + + if (dim == 2) + { + if ((coord[0] >= 0) + && (coord[0] < size[0]) + && (coord[1] >= 0 ) + && (coord[1] < size[1] )) + { + return true; + } + } + if (dim == 3) + { + if ((coord[0] >= 0) + && (coord[0] < size[0]) + && (coord[1] >= 0 ) + && (coord[1] < size[1] ) + && (coord[2] >= 0 ) + && (coord[2] < size[2] )) + { + return true; + } + } + return false; + } + + + template + inline std::vector< ShortestPathNode* > + ShortestPathImageFilter:: + GetNeighbors (unsigned int nodeNum, bool FullNeighbors) + { + // returns a vector of nodepointers.. these nodes are the neighbors + int dim = InputImageType::ImageDimension; + IndexType Coord = NodeToCoord(nodeNum); + IndexType NeighborCoord; + std::vector nodeList; + + int neighborDistance = 1; //if i increase that, i might not hit the endnote + + // maybe use itkNeighborhoodIterator here, might be faster + + if ( dim == 2) + { + // N4 + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + if (FullNeighbors) + { + // N8 + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + } + } + if ( dim == 3) + { + // N6 + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + if (FullNeighbors) + { + // N26 + // Middle Slice + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + // BackSlice (Diagonal) + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + //BackSlice (Non-Diag) + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]-neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + // FrontSlice (Diagonal) + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + //FrontSlice(Non-Diag) + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]-neighborDistance; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]+neighborDistance; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]; + NeighborCoord[1] = Coord[1]+neighborDistance; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + NeighborCoord[0] = Coord[0]-neighborDistance; + NeighborCoord[1] = Coord[1]; + NeighborCoord[2] = Coord[2]+neighborDistance; + if (CoordIsInBounds(NeighborCoord)) + nodeList.push_back(&m_Nodes[CoordToNode(NeighborCoord)]); + + } + } + return nodeList; + } + + + template + void ShortestPathImageFilter:: + SetStartIndex (const typename TInputImageType::IndexType &StartIndex) + { + for (unsigned int i=0;i v; + v[0] = m_EndIndex[0]-a[0]; + v[1] = m_EndIndex[1]-a[1]; + v[2] = m_EndIndex[2]-a[2]; + + return m_CostFunction->GetMinCost() * v.GetNorm(); + } + + + template + void + ShortestPathImageFilter:: + InitGraph() + { + // Clean up previous stuff + CleanUp(); + + // initalize cost function + m_CostFunction->Initialize(); + + // Calc Number of nodes + m_ImageDimensions = TInputImageType::ImageDimension; + const InputImageSizeType &size = this->GetInput()->GetRequestedRegion().GetSize(); + m_Graph_NumberOfNodes = 1; + for (NodeNumType i=0; i + void + ShortestPathImageFilter:: + StartShortestPathSearch() + { + // Setup Timer + clock_t startAll = clock(); + clock_t stopAll = clock(); + + // init variables + double durationAll = 0; + bool timeout = false; + bool makeNewHeapNecessary = false; + NodeNumType mainNodeListIndex = 0; + DistanceType curNodeDistance = 0; + DistanceType curNodeDistAndEst = 0; + NodeNumType numberOfNodesChecked = 0; + + // Create Multimap (tree structure for fast searching) + std::multimap myMap; + std::pair< std::multimap::iterator, std::multimap::iterator> ret; + std::multimap::iterator it; + + // At first, only startNote is discovered. + myMap.insert( std::pair (m_Nodes[m_Graph_StartNode].distAndEst, &m_Nodes[m_Graph_StartNode]) ); + + // While there are discovered Nodes, pick the one with lowest distance, + // update its neighbors and eventually delete it from the discovered Nodes list. + while(!myMap.empty()) + { + numberOfNodesChecked++; + if ( (numberOfNodesChecked % (m_Graph_NumberOfNodes/100)) == 0) + { + MITK_INFO << "Checked " << ( numberOfNodesChecked / (m_Graph_NumberOfNodes/100) ) << "% List: " << myMap.size() << "\n"; + } + + // Get element with lowest score + mainNodeListIndex = myMap.begin()->second->mainListIndex; + curNodeDistAndEst = myMap.begin()->second->distAndEst; + curNodeDistance = myMap.begin()->second->distance; + myMap.begin()->second->closed = true; // close it + + // Debug: + //MITK_INFO << "INFO: size " << myMap.size(); + /* + for (it = myMap.begin(); it != myMap.end(); ++it) + { + MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex; + } + */ + + // Kicks out element with lowest score + myMap.erase( myMap.begin() ); + + // if wanted, store vector order + if (m_StoreVectorOrder) + { + m_VectorOrder.push_back(mainNodeListIndex); + } + + // Check neighbors + std::vector neighborNodes = GetNeighbors(mainNodeListIndex, m_Graph_fullNeighbors); + for (NodeNumType i=0; iclosed) + continue; // this nodes is already closed, go to next neighbor + + IndexType coordCurNode = NodeToCoord(mainNodeListIndex); + IndexType coordNeighborNode = NodeToCoord(neighborNodes[i]->mainListIndex); + + // calculate the new Distance to the current neighbor + double newDistance = curNodeDistance + + (m_CostFunction->GetCost(coordCurNode, coordNeighborNode)); + + // if it is shorter than any yet known path to this neighbor, than the current path is better. Save that! + if ((newDistance < neighborNodes[i]->distance) || (neighborNodes[i]->distance == -1) ) + { + // if that neighbornode is not in discoverednodeList yet, Push it there and update + if (neighborNodes[i]->distance == -1) + { + + neighborNodes[i]->distance = newDistance; + neighborNodes[i]->distAndEst = newDistance + getEstimatedCostsToTarget(coordNeighborNode); + neighborNodes[i]->prevNode = mainNodeListIndex; + myMap.insert( std::pair (m_Nodes[neighborNodes[i]->mainListIndex].distAndEst, &m_Nodes[neighborNodes[i]->mainListIndex]) ); + /* + MITK_INFO << "Inserted: " << m_Nodes[neighborNodes[i]->mainListIndex].distAndEst << "|" << m_Nodes[neighborNodes[i]->mainListIndex].mainListIndex; + MITK_INFO << "INFO: size " << myMap.size(); + for (it = myMap.begin(); it != myMap.end(); ++it) + { + MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex; + } + */ + } + // or if is already in discoverednodelist, update + else + { + + /* + it = myMap.find(neighborNodes[i]->distAndEst); + if (it == myMap.end() ) + { + MITK_INFO << "Nothing!"; + // look further + for (it = myMap.begin(); it != myMap.end(); ++it) + { + if ((*it).second->mainListIndex == lookForId) + { + MITK_INFO << "But it is there!!!"; + MITK_INFO << "Searched for: " << lookFor << " but had: " << (*it).second->distAndEst; + } + + } + } + */ + + // 1st : find and delete old element + bool found = false; + double lookFor = neighborNodes[i]->distAndEst; + unsigned int lookForId = neighborNodes[i]->mainListIndex; + ret = myMap.equal_range(neighborNodes[i]->distAndEst); + + if ((ret.first == ret.second)) + { + MITK_INFO << "No exact match!"; // if this happens, you are screwed + /* + MITK_INFO << "Was looking for: " << lookFor << " ID: " << lookForId; + if (ret.first != myMap.end() ) + { + it = ret.first; + MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex; + ++it; + MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex; + --it; + --it; + MITK_INFO << "Found: " << it->first << " ID: " << it->second->mainListIndex; + } + + // look if that ID is found in the map + for (it = myMap.begin(); it != myMap.end(); ++it) + { + if ((*it).second->mainListIndex == lookForId) + { + MITK_INFO << "But it is there!!!"; + MITK_INFO << "Searched dist: " << lookFor << " found dist: " << (*it).second->distAndEst; + } + + } + */ + } + else + { + for (it=ret.first; it!=ret.second; ++it) + { + if (it->second->mainListIndex == neighborNodes[i]->mainListIndex) + { + found = true; + myMap.erase(it); + /* + MITK_INFO << "INFO: size " << myMap.size(); + MITK_INFO << "Erase: " << it->first << "|" << it->second->mainListIndex; + + MITK_INFO << "INFO: size " << myMap.size(); + for (it = myMap.begin(); it != myMap.end(); ++it) + { + MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex; + } + */ + break; + } + } + } + + + if (!found) + { + MITK_INFO << "Could not find it! :("; + continue; + } + + // 2nd: update and insert new element + neighborNodes[i]->distance = newDistance; + neighborNodes[i]->distAndEst = newDistance + getEstimatedCostsToTarget(coordNeighborNode); + neighborNodes[i]->prevNode = mainNodeListIndex; + //myMap.insert( std::pair (neighborNodes[i]->distAndEst, neighborNodes[i])); + myMap.insert( std::pair (m_Nodes[neighborNodes[i]->mainListIndex].distAndEst, &m_Nodes[neighborNodes[i]->mainListIndex]) ); + + //MITK_INFO << "Re-Inserted: " << m_Nodes[neighborNodes[i]->mainListIndex].distAndEst << "|" << m_Nodes[neighborNodes[i]->mainListIndex].mainListIndex; + //MITK_INFO << "INFO: size " << myMap.size(); + /*for (it = myMap.begin(); it != myMap.end(); ++it) + { + MITK_INFO << "(1) " << it->first << "|" << it->second->distAndEst << "|"<second->mainListIndex; + }*/ + + } + } + } + // finished with checking all neighbors. + + // Check Timeout, if activated + if (m_ActivateTimeOut) + { + stopAll = clock(); + durationAll = (double)(stopAll - startAll) / CLOCKS_PER_SEC; + if (durationAll >= 30) + { + MITK_INFO << "TIMEOUT!! Search took over 30 seconds"; + timeout = true ; + } + } + + // Check end criteria: + // For multiple points + if ( multipleEndPoints ) + { + // super slow, make it faster + for (int i=0 ;i + void + ShortestPathImageFilter:: + MakeOutputs() + { + MITK_INFO << "Make Output"; + + if (m_MakeOutputImage) + { + OutputImagePointer output0 = this->GetOutput(0); + output0->SetRegions( this->GetInput()->GetLargestPossibleRegion() ); + output0->Allocate(); + OutputImageIteratorType shortestPathImageIt (output0, output0->GetRequestedRegion()); + // Create ShortestPathImage (Output 0) + for (shortestPathImageIt.GoToBegin(); !shortestPathImageIt.IsAtEnd(); ++shortestPathImageIt) + { + // First intialize with background color + shortestPathImageIt.Set( BACKGROUND ) ; + } + + if (!multipleEndPoints) // Only one path was calculated + { + for (int i=0; i< m_VectorPath.size(); i++) + { + shortestPathImageIt.SetIndex( m_VectorPath[i] ); + shortestPathImageIt.Set( FOREGROUND ) ; + } + } + else // multiple pathes has been calculated, draw all + { + for (int i =0; i + typename ShortestPathImageFilter::OutputImagePointer + ShortestPathImageFilter:: + GetVectorOrderImage() + { + // Create Vector Order Image + // Return it + + OutputImagePointer image = OutputImageType::New(); + image->SetRegions( this->GetInput()->GetLargestPossibleRegion() ); + image->Allocate(); + OutputImageIteratorType vectorOrderImageIt (image, image->GetRequestedRegion()); + + + MITK_INFO << "GetVectorOrderImage"; + for (vectorOrderImageIt.GoToBegin(); !vectorOrderImageIt.IsAtEnd(); ++vectorOrderImageIt) + { + // First intialize with background color + vectorOrderImageIt.Value() = BACKGROUND ; + } + for (int i=0; i< m_VectorOrder.size(); i++) + { + vectorOrderImageIt.SetIndex(NodeToCoord(m_VectorOrder[i]) ); + vectorOrderImageIt.Set( BACKGROUND+i) ; + } + return image; + } + + + template + typename ShortestPathImageFilter::OutputImagePointer + ShortestPathImageFilter:: + GetDistanceImage() + { + // Create Distance Image + // Return it + + OutputImagePointer image = OutputImageType::New(); + image->SetRegions( this->GetInput()->GetLargestPossibleRegion() ); + image->Allocate();; + OutputImageIteratorType distanceImageIt (image, image->GetRequestedRegion()); + // Create Distance Image (Output 1) + NodeNumType myNodeNum; + for (distanceImageIt.GoToBegin(); !distanceImageIt.IsAtEnd(); ++distanceImageIt) + { + IndexType index = distanceImageIt.GetIndex(); + myNodeNum = CoordToNode(index); + double newVal = m_Nodes[myNodeNum].distance; + distanceImageIt.Set(newVal); + } + } + + + + template + std::vector< itk::Index<3> > + ShortestPathImageFilter:: + GetVectorPath() + { + return m_VectorPath; + } + + template + std::vector< std::vector< itk::Index<3> > > + ShortestPathImageFilter:: + GetMultipleVectorPaths() + { + return m_MultipleVectorPaths; + } + + template + void + ShortestPathImageFilter:: + MakeShortestPathVector() + { + MITK_INFO << "Make ShortestPath Vec"; + + // single end point + if ( !multipleEndPoints ) + { + // fill m_VectorPath with the Shortest Path + m_VectorPath.clear(); + + // Go backwards from endnote to startnode + NodeNumType prevNode = m_Graph_EndNode; + while (prevNode != -1) + { + m_VectorPath.push_back( NodeToCoord(prevNode) ); + + if (prevNode == m_Graph_StartNode) + break; + + prevNode = m_Nodes[prevNode].prevNode; + } + + // reverse it + std::reverse(m_VectorPath.begin(), m_VectorPath.end() ); + } + // Multiple end end points and pathes + else + { + for (int i=0; i + void + ShortestPathImageFilter:: + CleanUp() + { + + + m_VectorPath.clear(); + //TODO: if multiple Path, clear all multiple Paths + /* + for (NodeNumType i=0; i + void + ShortestPathImageFilter:: + GenerateData() + { + // Build Graph + InitGraph(); + + // Calc Shortest Parth + StartShortestPathSearch(); + + // Fill Shortest Path + MakeShortestPathVector(); + + // Make Outputs + MakeOutputs(); + } + + template + void + ShortestPathImageFilter:: + PrintSelf( std::ostream& os, Indent indent ) const + { + Superclass::PrintSelf(os,indent); + + } + +} /* end namespace itk */ diff --git a/Modules/GraphAlgorithms/itkShortestPathNode.cpp b/Modules/GraphAlgorithms/itkShortestPathNode.cpp new file mode 100644 index 0000000000..975a975885 --- /dev/null +++ b/Modules/GraphAlgorithms/itkShortestPathNode.cpp @@ -0,0 +1,15 @@ +#include "itkShortestPathNode.h" + +namespace itk +{ + +// bool ShortestPathNode::operator<(const ShortestPathNode &a) const +// { +// return (this->distance < a.distance); +// } +// bool ShortestPathNode::operator==(const ShortestPathNode &a) const +// { +// return (this->mainListIndex == a.mainListIndex); +// } + +} diff --git a/Modules/GraphAlgorithms/itkShortestPathNode.h b/Modules/GraphAlgorithms/itkShortestPathNode.h new file mode 100644 index 0000000000..787b03f0a1 --- /dev/null +++ b/Modules/GraphAlgorithms/itkShortestPathNode.h @@ -0,0 +1,27 @@ +#ifndef __itkShortestPathNode_h_ +#define __itkShortestPathNode_h_ + + +namespace itk +{ + typedef double DistanceType; // Type to declare the costs + typedef unsigned int NodeNumType; // Type for Node Numeration: unsignend int for up to 4.2 billion pixel in 32Bit system + + class ShortestPathNode + { + public: + DistanceType distance; // minimal costs from StartPoint to this pixel + DistanceType distAndEst; // Distance+Estimated Distnace to target + NodeNumType prevNode; // previous node. Important to find the Shortest Path + NodeNumType mainListIndex; // Indexnumber of this node in m_Nodes + bool closed; // determines if this node is closes, so its optimal path to startNode is known + }; + + //bool operator<(const ShortestPathNode &a) const; + //bool operator==(const ShortestPathNode &a) const; + +} + + + +#endif