diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp index bfa25499af..7c408361cc 100644 --- a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp +++ b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp @@ -1,212 +1,209 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkTractAnalyzer_cpp #define __mitkTractAnalyzer_cpp #include #include #include #include #include #include #include #include #include using namespace std; namespace mitk { - TractAnalyzer::TractAnalyzer() - { - - } + TractAnalyzer::TractAnalyzer() { } - void TractAnalyzer::BuildGraph() + void TractAnalyzer::MakeRoi() { 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; iProgress(); mitk::Point3D p = m_PointSetNode->GetPoint(i); mitk::Point3D p2 = m_PointSetNode->GetPoint(i+1); itk::Index<3> startPoint; itk::Index<3> endPoint; m_InputImage->GetGeometry()->WorldToIndex(p,startPoint); m_InputImage->GetGeometry()->WorldToIndex(p2,endPoint); MITK_INFO << "create roi"; std::vector< itk::Index<3> > path = CreateSegment(startPoint, endPoint); for(std::vector< itk::Index<3> >::iterator it = path.begin(); it != path.end(); it++) { itk::Index<3> ix = *it; if (!(ix==endPoint)) { mitk::ProgressBar::GetInstance()->Progress(); 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(); } } } } // save pathDescription to m_PathDescription m_PathDescription = pathDescription; FloatImageType::Pointer itkImg = FloatImageType::New(); mitk::CastToItkImage(m_InputImage, 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); } m_TbssRoi = mitk::TbssRoiImage::New(); m_TbssRoi->SetRoi(roi); m_TbssRoi->SetImage(roiImg); m_TbssRoi->InitializeFromImage(); } } std::vector< itk::Index<3> > TractAnalyzer::CreateSegment(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(); return pathFinder->GetVectorPath(); } } } #endif diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h index 09e543dc6b..f1f0041c1e 100644 --- a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h +++ b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h @@ -1,125 +1,148 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkTractAnalyzer_h_ #define __mitkTractAnalyzer_h_ #include "QuantificationExports.h" #include -#include #include "mitkImage.h" #include "mitkImageCast.h" #include #include namespace mitk{ /** \class TractAnalyzer */ class Quantification_EXPORT TractAnalyzer { public: TractAnalyzer(); - ~TractAnalyzer() {}; + ~TractAnalyzer() {} + + /** Image type definitions */ typedef itk::Image CharImageType; typedef itk::Image FloatImageType; - typedef itk::Image ProjectionsImageType; - typedef itk::VectorImage VectorImageType; - - void BuildGraph(); - - - /* - std::vector< itk::Index<3> > GetPath() - { - return m_Path; - } - - CharImageType::Pointer GetRoiImage() - { - return m_RoiImg; - } + /** \brief Main method for region of interest calculation + * + * A region of interest is calculated adding the segments between the points on the ROI + * that was specified by the user. + */ + void MakeRoi(); -*/ + /** \brief Sets the input image + * + * The region of interest is calculated on a 3D image. This is generally the mean FA skeleton as calculated + * in the standard TBSS pipeline (see http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TBSS). + */ mitk::TbssRoiImage::Pointer GetRoiImage() { return m_TbssRoi; } + /** \brief Sets the input image + * + * The region of interest is calculated on a 3D image. This is generally the mean FA skeleton as calculated + * in the standard TBSS pipeline (see http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TBSS). + */ void SetInputImage(mitk::Image::Pointer inputImage) { m_InputImage = inputImage; } + + /** \brief Sets the user-defined point set + * + * Set the user-defined point sets. The region of interest must pass through these points. + */ void SetPointSet(mitk::PointSet::Pointer pointSet) { m_PointSetNode = pointSet; } + /** \brief Sets a lower bound for the threshold. + * + * Low fractional anisotropy values can indicate partial volume of non white matter tissue. + * This thresholds limits the search for a region of interest to voxels with a minimum value. + */ void SetThreshold(double threshold) { m_Threshold = threshold; } + + + /** \brief Returns a string with the indices of points on the region of interest + * + * The region of interest calculated by the TractAnalyzer contains a list of ITK indices. + * This method returns a string containing these indices for display in the GUI + */ std::string GetPathDescription() { return m_PathDescription; } protected: + /** \brief Calculates a segment of the region of interest + * + * The region of interest is calculated on a 3D image. This is generally the mean FA skeleton as calculated + * in the standard TBSS pipeline (see http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TBSS). + */ std::vector< itk::Index<3> > CreateSegment(itk::Index<3> startPoint, itk::Index<3> endPoint); + /** \brief Output TbssRoiImage */ mitk::TbssRoiImage::Pointer m_TbssRoi; - //CharImageType::Pointer m_RoiImg; - + /** \brief Inputimage */ mitk::Image::Pointer m_InputImage; + /** \brief Threshold for ROI search */ double m_Threshold; - //std::vector< itk::Index<3> > m_Path; - + /** \brief User defined point set */ mitk::PointSet::Pointer m_PointSetNode; + /** \brief Path description in as string for display in GUI */ std::string m_PathDescription; private: }; } #endif //__itkTractAnalyzer_h_ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp index 79bc0e16dc..1bcbc2b222 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp @@ -1,1315 +1,1315 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // 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 "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 #include #include "mitkITKImageImport.h" // #include "mitkImageMapperGL2D.h" #include "mitkVolumeDataVtkMapper3D.h" #include "mitkImageAccessByItk.h" #include "mitkTensorImage.h" #include "itkDiffusionTensor3D.h" #define SEARCHSIGMA 10 // length in linear voxel dimens #define MAXSEARCHLENGTH (3*SEARCHSIGMA) const std::string QmitkTractbasedSpatialStatisticsView::VIEW_ID = "org.mitk.views.tractbasedspatialstatistics"; using namespace berry; QmitkTractbasedSpatialStatisticsView::QmitkTractbasedSpatialStatisticsView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) { } QmitkTractbasedSpatialStatisticsView::~QmitkTractbasedSpatialStatisticsView() { } void QmitkTractbasedSpatialStatisticsView::PerformChange() { m_Controls->m_RoiPlotWidget->ModifyPlot(m_Controls->m_Segments->value(), m_Controls->m_Average->isChecked()); } 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"); } bool foundTbssRoi = false; bool foundTbss = false; bool found3dImage = false; bool found4dImage = false; bool foundFiberBundle = false; bool foundStartRoi = false; bool foundEndRoi = false; mitk::TbssRoiImage* roiImage; mitk::TbssImage* image; mitk::Image* img; mitk::FiberBundleX* fib; mitk::PlanarFigure* start; mitk::PlanarFigure* end; m_CurrentStartRoi = NULL; m_CurrentEndRoi = NULL; for ( int i=0; iGetData(); if( nodeData ) { if(QString("TbssRoiImage").compare(nodeData->GetNameOfClass())==0) { foundTbssRoi = true; roiImage = static_cast(nodeData); } else if (QString("TbssImage").compare(nodeData->GetNameOfClass())==0) { foundTbss = true; image = static_cast(nodeData); } else if(QString("Image").compare(nodeData->GetNameOfClass())==0) { img = static_cast(nodeData); if(img->GetDimension() == 3) { found3dImage = true; } else if(img->GetDimension() == 4) { found4dImage = true; } } else if (QString("FiberBundleX").compare(nodeData->GetNameOfClass())==0) { foundFiberBundle = true; fib = static_cast(nodeData); this->m_CurrentFiberNode = nodes[i]; } if(QString("PlanarCircle").compare(nodeData->GetNameOfClass())==0) { if(!foundStartRoi) { start = dynamic_cast(nodeData); this->m_CurrentStartRoi = nodes[i]; foundStartRoi = true; } else { end = dynamic_cast(nodeData); this->m_CurrentEndRoi = nodes[i]; foundEndRoi = true; } } } } this->m_Controls->m_CreateRoi->setEnabled(found3dImage); this->m_Controls->m_ImportFsl->setEnabled(found4dImage); if(foundTbss && foundTbssRoi) { this->Plot(image, roiImage); } if(found3dImage && foundFiberBundle && foundStartRoi && foundEndRoi) { this->PlotFiberBundle(fib, img, start, end); } else if(found3dImage == true && foundFiberBundle) { this->PlotFiberBundle(fib, img); } if(found3dImage) { this->InitPointsets(); } this->m_Controls->m_Cut->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi); this->m_Controls->m_SegmentLabel->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); this->m_Controls->m_Segments->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); this->m_Controls->m_Average->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); } 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(true) ); // 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_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&) ) ); connect( (QObject*)(m_Controls->m_Cut), SIGNAL(clicked()), this, SLOT(Cut()) ); connect( (QObject*)(m_Controls->m_Average), SIGNAL(stateChanged(int)), this, SLOT(PerformChange()) ); connect( (QObject*)(m_Controls->m_Segments), SIGNAL(valueChanged(int)), this, SLOT(PerformChange()) ); } } void QmitkTractbasedSpatialStatisticsView::CopyToClipboard() { if(m_Controls->m_RoiPlotWidget->IsPlottingFiber()) { // Working with fiber bundles std::vector > profiles = m_Controls->m_RoiPlotWidget->GetIndividualProfiles(); QString clipboardText; for (std::vector >::iterator it = profiles.begin(); it != profiles.end(); ++it) { for (std::vector::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); } clipboardText.append(QString("\n")); } if(m_Controls->m_Average->isChecked()) { std::vector averages = m_Controls->m_RoiPlotWidget->GetAverageProfile(); clipboardText.append(QString("\nAverage\n")); for (std::vector::iterator it2 = averages.begin(); it2 != averages.end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); } } QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); } else{ // Working with TBSS Data if(m_Controls->m_Average->isChecked()) { 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); } else { std::vector > vals = m_Controls->m_RoiPlotWidget->GetIndividualProfiles(); 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 = ""; std::vector nodes = this->GetDataManagerSelection(); for ( int i=0; iGetData()->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodes[i]->GetData()); if(img->GetDimension() == 4) { importer->SetImportVolume(img); name = nodes[i]->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) { int index = (int)pos.x(); if(m_Roi.size() > 0 && m_CurrentGeometry != NULL && !m_Controls->m_RoiPlotWidget->IsPlottingFiber() ) { index = std::min( (int)m_Roi.size()-1, std::max(0, index) ); 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); } else if(m_Controls->m_RoiPlotWidget->IsPlottingFiber() ) { mitk::Point3D point = m_Controls->m_RoiPlotWidget->GetPositionInWorld(index); m_MultiWidget->MoveCrossToPosition(point); } } void QmitkTractbasedSpatialStatisticsView::Cut() { mitk::BaseData* fibData = m_CurrentFiberNode->GetData(); mitk::FiberBundleX* fib = static_cast(fibData); mitk::BaseData* startData = m_CurrentStartRoi->GetData(); mitk::PlanarFigure* startRoi = static_cast(startData); mitk::PlaneGeometry* startGeometry2D = dynamic_cast( const_cast(startRoi->GetGeometry2D()) ); mitk::BaseData* endData = m_CurrentEndRoi->GetData(); mitk::PlanarFigure* endRoi = static_cast(endData); mitk::PlaneGeometry* endGeometry2D = dynamic_cast( const_cast(endRoi->GetGeometry2D()) ); mitk::Point3D startCenter = startRoi->GetWorldControlPoint(0); //center Point of start roi mitk::Point3D endCenter = endRoi->GetWorldControlPoint(0); //center Point of end roi mitk::FiberBundleX::Pointer inStart = fib->ExtractFiberSubset(startRoi); mitk::FiberBundleX::Pointer inBoth = inStart->ExtractFiberSubset(endRoi); int num = inBoth->GetNumFibers(); vtkSmartPointer fiberPolyData = inBoth->GetFiberPolyData(); vtkCellArray* lines = fiberPolyData->GetLines(); lines->InitTraversal(); // initialize new vtk polydata vtkSmartPointer points = vtkSmartPointer::New(); vtkSmartPointer polyData = vtkSmartPointer::New(); vtkSmartPointer cells = vtkSmartPointer::New(); int pointIndex=0; // find start and endpoint for( int fiberID( 0 ); fiberID < num; fiberID++ ) { vtkIdType numPointsInCell(0); vtkIdType* pointsInCell(NULL); lines->GetNextCell ( numPointsInCell, pointsInCell ); int startId = 0; int endId = 0; float minDistStart = std::numeric_limits::max(); float minDistEnd = std::numeric_limits::max(); vtkSmartPointer polyLine = vtkSmartPointer::New(); int lineIndex=0; for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); mitk::Point3D point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; float distanceToStart = point.EuclideanDistanceTo(startCenter); float distanceToEnd = point.EuclideanDistanceTo(endCenter); if(distanceToStart < minDistStart) { minDistStart = distanceToStart; startId = pointInCellID; } if(distanceToEnd < minDistEnd) { minDistEnd = distanceToEnd; endId = pointInCellID; } } /* We found the start and end points of of the part that should be plottet for the current fiber. now we need to plot them. If the endId is smaller than the startId the plot order must be reversed*/ if(startId < endId) { double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId+1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool secondOnPositive = startGeometry2D->IsAbove(pSecond); // Calculate intersection with the plane mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId+1 lie on the same side of the plane, so we need need startId-1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; double insertPoint[3]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; // First insert the intersection with the start roi points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; if(! (secondOnPositive ^ startOnPositive) ) { /* StartId and startId+1 lie on the same side of the plane so startId is also part of the ROI*/ double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); points->InsertNextPoint(start); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Insert the rest up and to including endId-1 for( int pointInCellID( startId+1 ); pointInCellID < endId ; pointInCellID++) { // create new polyline for new polydata double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); points->InsertNextPoint(p); // add point to line polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } /* endId must be included if endId and endId-1 lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; bool lastOnPositive = endGeometry2D->IsAbove(pLast); bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); normal = endGeometry2D->GetNormal(); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId-1 lie on the same side of the plane, so we need need endId+1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure endId is part of the roi so can also be included here*/ p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; double *end = fiberPolyData->GetPoint( pointsInCell[endId] ); points->InsertNextPoint(end); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; //Insert the Last Point (intersection with the end roi) points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Need to reverse walking order else{ double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool secondOnPositive = startGeometry2D->IsAbove(pSecond); // Calculate intersection with the plane mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId-1 lie on the same side of the plane, so we need need startId+1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; double insertPoint[3]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; // First insert the intersection with the start roi points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId-1 lie on the same side of the plane so endId is also part of the ROI*/ double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); points->InsertNextPoint(start); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Insert the rest up and to including endId-1 for( int pointInCellID( startId-1 ); pointInCellID > endId ; pointInCellID--) { // create new polyline for new polydata double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); points->InsertNextPoint(p); // add point to line polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } /* startId must be included if startId and startId+ lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; bool lastOnPositive = endGeometry2D->IsAbove(pLast); mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId+1 lie on the same side of the plane, so we need need endId-1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; /* endId and endId+1 lie on the same side of the plane so startId is also part of the ROI*/ double *end = fiberPolyData->GetPoint( pointsInCell[endId] ); points->InsertNextPoint(end); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; //Insert the Last Point (intersection with the end roi) points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // add polyline to vtk cell array cells->InsertNextCell(polyLine); } // Add the points to the dataset polyData->SetPoints(points); // Add the lines to the dataset polyData->SetLines(cells); mitk::FiberBundleX::Pointer cutBundle = mitk::FiberBundleX::New(polyData); mitk::DataNode::Pointer cutNode = mitk::DataNode::New(); cutNode->SetData(cutBundle); std::string name = "fiberCut"; cutNode->SetName(name); GetDataStorage()->Add(cutNode); } 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::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 bool ok; double threshold = QInputDialog::getDouble(m_Controls->m_CreateRoi, tr("Set an FA threshold"), tr("Threshold:"), 0.2, 0.0, 1.0, 2, &ok); if(!ok) return; mitk::Image::Pointer image; std::vector nodes = this->GetDataManagerSelection(); for ( int i=0; iGetData()->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodes[i]->GetData()); if(img->GetDimension() == 3) { image = img; } } } if(image.IsNull()) { return; } mitk::TractAnalyzer analyzer; analyzer.SetInputImage(image); analyzer.SetThreshold(threshold); // Set Pointset to analyzer analyzer.SetPointSet(m_PointSetNode); // Run Analyzer - analyzer.BuildGraph(); + analyzer.MakeRoi(); // Obtain tbss roi image from analyzer mitk::TbssRoiImage::Pointer tbssRoi = analyzer.GetRoiImage(); tbssRoi->SetStructure(m_Controls->m_Structure->text().toStdString()); // get path description and set to interface std::string pathDescription = analyzer.GetPathDescription(); m_Controls->m_PathTextEdit->setPlainText(QString(pathDescription.c_str())); // Add roi image to datastorage AddTbssToDataStorage(tbssRoi, m_Controls->m_RoiName->text().toStdString()); } void QmitkTractbasedSpatialStatisticsView:: PlotFiberBundle(mitk::FiberBundleX *fib, mitk::Image* img, mitk::PlanarFigure* startRoi, mitk::PlanarFigure* endRoi) { bool avg = m_Controls->m_Average->isChecked(); int segments = m_Controls->m_Segments->value(); m_Controls->m_RoiPlotWidget->PlotFiberBetweenRois(fib, img, startRoi ,endRoi, avg, segments); m_Controls->m_RoiPlotWidget->SetPlottingFiber(true); mitk::RenderingManager::GetInstance()->ForceImmediateUpdateAll(); } 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 = ""; std::string structure = roiImage->GetStructure(); m_Controls->m_RoiPlotWidget->SetGroups(image->GetGroupInfo()); 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); } m_Controls->m_RoiPlotWidget->SetPlottingFiber(false); } 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; } }