diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp index e7220cdb35..8d48e6f0f7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/QmitkTbssRoiAnalysisWidget.cpp @@ -1,967 +1,967 @@ /*=================================================================== 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. ===================================================================*/ #include "QmitkTbssRoiAnalysisWidget.h" #include "mitkImagePixelReadAccessor.h" #include "mitkPixelTypeMultiplex.h" #include #include #include #include QmitkTbssRoiAnalysisWidget::QmitkTbssRoiAnalysisWidget( QWidget * parent ) : QmitkPlotWidget(parent) { m_PlotPicker = new QwtPlotPicker(m_Plot->canvas()); m_PlotPicker->setStateMachine(new QwtPickerDragPointMachine()); m_PlotPicker->setTrackerMode(QwtPicker::ActiveOnly); m_PlottingFiberBundle = false; } void QmitkTbssRoiAnalysisWidget::DoPlotFiberBundles(mitk::FiberBundle *fib, mitk::Image* img, mitk::DataNode* startRoi, mitk::DataNode* endRoi, bool avg, int number) { TractContainerType tracts = CreateTracts(fib, startRoi, endRoi); TractContainerType resampledTracts = ParameterizeTracts(tracts, number); // Now we have the resampled tracts. Next we should use these points to read out the values mitkPixelTypeMultiplex3(PlotFiberBundles,img->GetImageDescriptor()->GetChannelTypeById(0),resampledTracts, img, avg); m_CurrentTracts = resampledTracts; } TractContainerType QmitkTbssRoiAnalysisWidget::CreateTracts(mitk::FiberBundle *fib, mitk::DataNode *startRoi, mitk::DataNode *endRoi) { mitk::PlaneGeometry* startGeometry2D = const_cast(dynamic_cast(startRoi->GetData())->GetPlaneGeometry()); mitk::PlaneGeometry* endGeometry2D = const_cast(dynamic_cast(endRoi->GetData())->GetPlaneGeometry()); mitk::Point3D startCenter = dynamic_cast(startRoi->GetData())->GetWorldControlPoint(0); //center Point of start roi mitk::Point3D endCenter = dynamic_cast(startRoi->GetData())->GetWorldControlPoint(0); //center Point of end roi mitk::FiberBundle::Pointer inStart = fib->ExtractFiberSubset(startRoi, nullptr); mitk::FiberBundle::Pointer inBoth = inStart->ExtractFiberSubset(endRoi, nullptr); int num = inBoth->GetNumFibers(); TractContainerType tracts; vtkSmartPointer fiberPolyData = inBoth->GetFiberPolyData(); vtkCellArray* lines = fiberPolyData->GetLines(); lines->InitTraversal(); // Now find out for each fiber which ROI is encountered first. If this is the startRoi, the direction is ok // Otherwise the plot should be in the reverse direction for( int fiberID( 0 ); fiberID < num; fiberID++ ) { vtkIdType numPointsInCell(0); vtkIdType* pointsInCell(nullptr); lines->GetNextCell ( numPointsInCell, pointsInCell ); int startId = 0; int endId = 0; mitk::ScalarType minDistStart = std::numeric_limits::max(); mitk::ScalarType minDistEnd = std::numeric_limits::max(); for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); mitk::Point3D point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; mitk::ScalarType distanceToStart = point.EuclideanDistanceTo(startCenter); mitk::ScalarType 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*/ TractType singleTract; PointType point; if(startId < endId) { // Calculate the intersection of the ROI with the startRoi and decide if the startId is part of the roi or must be cut of mitk::ScalarType *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]; mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); bool secondOnPositive = startGeometry2D->IsAbove(pSecond); 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]; } mitk::ScalarType d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); if(! (secondOnPositive ^ startOnPositive) ) { /* StartId and startId+1 lie on the same side of the plane so startId is also part of the ROI*/ mitk::ScalarType *start = fiberPolyData->GetPoint( pointsInCell[startId] ); point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; singleTract.push_back(point); } for( int pointInCellID( startId+1 ); pointInCellID < endId ; pointInCellID++) { // push back point mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); } /* 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]; normal = endGeometry2D->GetNormal(); bool lastOnPositive = endGeometry2D->IsAbove(pLast); 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 endId is also part of the tract and can be inserted here */ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); } else{ // Calculate the intersection of the ROI with the startRoi and decide if the startId is part of the roi or must be cut of mitk::ScalarType *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]; mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); bool secondOnPositive = startGeometry2D->IsAbove(pSecond); 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]; } mitk::ScalarType d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); if(! (secondOnPositive ^ startOnPositive) ) { /* StartId and startId+1 lie on the same side of the plane so startId is also part of the ROI*/ mitk::ScalarType *start = fiberPolyData->GetPoint( pointsInCell[startId] ); point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; singleTract.push_back(point); } for( int pointInCellID( startId-1 ); pointInCellID > endId ; pointInCellID--) { // push back point mitk::ScalarType *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); } /* endId must be included if endId and endI+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 also part of the tract and can be inserted here */ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; singleTract.push_back( point ); p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); point[0] = d*newPoint[0] + p0[0]; point[1] = d*newPoint[1] + p0[1]; point[2] = d*newPoint[2] + p0[2]; singleTract.push_back(point); } tracts.push_back(singleTract); } return tracts; } void QmitkTbssRoiAnalysisWidget::PlotFiberBetweenRois(mitk::FiberBundle *fib, mitk::Image* img, mitk::DataNode* startRoi, mitk::DataNode* endRoi, bool avg, int number) { if(fib == nullptr || img == nullptr || startRoi == nullptr || endRoi == nullptr) return; m_Fib = fib; m_CurrentImage = img; m_CurrentStartRoi = startRoi; m_CurrentEndRoi = endRoi; DoPlotFiberBundles(fib, img, startRoi, endRoi, avg, number); } void QmitkTbssRoiAnalysisWidget::ModifyPlot(int number, bool avg) { if(m_Fib == nullptr || m_CurrentTbssImage == nullptr || m_CurrentStartRoi == nullptr || m_CurrentEndRoi == nullptr) return; if(m_PlottingFiberBundle) { DoPlotFiberBundles(m_Fib, m_CurrentImage, m_CurrentStartRoi, m_CurrentEndRoi, avg, number); } else { PlotFiber4D(m_CurrentTbssImage, m_Fib, m_CurrentStartRoi, m_CurrentEndRoi, number); } } TractContainerType QmitkTbssRoiAnalysisWidget::ParameterizeTracts(TractContainerType tracts, int number) { TractContainerType resampledTracts; for(auto it = tracts.begin(); it != tracts.end(); ++it) { TractType resampledTract; TractType tract = *it; // Calculate the total length mitk::ScalarType totalLength = 0; if(tract.size() < 2) continue; PointType p0 = tract.at(0); for(unsigned int i = 1; i(resampledTract.size()) <= (number+1); position+=stepSize) { /* In case we walked to far we need to find the next segment we are on and on what relative position on that tract we are on. Small correction for rounding errors */ while(locationBetween > distance+0.001) { if(tractCounter == tract.size()) std::cout << "problem"; // Determine by what distance we are no on the next segment locationBetween = locationBetween - distance; p0 = p1; p1 = tract.at(tractCounter); tractCounter++; distance = p0.EuclideanDistanceTo(p1); } // Direction PointType::VectorType direction = p1-p0; direction.Normalize(); PointType newSample = p0 + direction*locationBetween; resampledTract.push_back(newSample); locationBetween += stepSize; } resampledTracts.push_back(resampledTract); } return resampledTracts; } mitk::Point3D QmitkTbssRoiAnalysisWidget::GetPositionInWorld(int index) { mitk::ScalarType xSum = 0.0; mitk::ScalarType ySum = 0.0; mitk::ScalarType zSum = 0.0; for(auto it = m_CurrentTracts.begin(); it!=m_CurrentTracts.end(); ++it) { TractType tract = *it; PointType p = tract.at(index); xSum += p[0]; ySum += p[1]; zSum += p[2]; } int number = m_CurrentTracts.size(); mitk::ScalarType xPos = xSum / number; mitk::ScalarType yPos = ySum / number; mitk::ScalarType zPos = zSum / number; mitk::Point3D pos; pos[0] = xPos; pos[1] = yPos; pos[2] = zPos; return pos; } std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfiles() { MITK_INFO << "make profiles!"; std::vector< std::vector > profiles; 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++; } profiles.push_back(profile); } m_IndividualProfiles = profiles; // 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 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(unsigned int i=0; i > groupProfiles = CalculateGroupProfiles(); Plot(groupProfiles); } void QmitkTbssRoiAnalysisWidget::Plot(std::vector > groupProfiles) { this->Clear(); m_Vals.clear(); std::vector v1; std::vector xAxis; for(unsigned 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++; } auto 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(); } std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfilesFibers(mitk::TbssImage::Pointer tbssImage, mitk::FiberBundle *fib, mitk::DataNode* startRoi, mitk::DataNode* endRoi, int number) { TractContainerType tracts = CreateTracts(fib, startRoi, endRoi); TractContainerType resampledTracts = ParameterizeTracts(tracts, number); int nTracts = resampledTracts.size(); this->Clear(); // For every group we have m fibers * n subjects of profiles to fill std::vector< std::vector > profiles; // calculate individual profiles by going through all n subjects int size = m_Projections->GetVectorLength(); for(int s=0; s profile; auto it = resampledTracts[t].begin(); while(it != resampledTracts[t].end()) { PointType p = *it; PointType index; tbssImage->GetGeometry()->WorldToIndex(p, index); itk::Index<3> ix; ix[0] = index[0]; ix[1] = index[1]; ix[2] = index[2]; // Get value from image profile.push_back(m_Projections->GetPixel(ix).GetElement(s)); it++; } profiles.push_back(profile); } } m_IndividualProfiles = profiles; // Now create the group averages (every group contains m fibers * n_i group members std::vector< std::pair >::iterator it; it = m_Groups.begin(); int c = 0; //the current profile number // Calculate the group averages std::vector< std::vector > groupProfiles; 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(unsigned int i=0; i > groupProfiles = CalculateGroupProfilesFibers(tbssImage, fib, startRoi, endRoi, number); Plot(groupProfiles); } template -void QmitkTbssRoiAnalysisWidget::PlotFiberBundles(const mitk::PixelType ptype, TractContainerType tracts, mitk::Image *img, bool avg) +void QmitkTbssRoiAnalysisWidget::PlotFiberBundles(const mitk::PixelType, TractContainerType tracts, mitk::Image *img, bool avg) { m_PlottingFiberBundle = true; this->Clear(); auto it = tracts.begin(); std::vector< std::vector > profiles; mitk::ImagePixelReadAccessor imAccess(img,img->GetVolumeData(0)); it = tracts.begin(); while(it != tracts.end()) { TractType tract = *it; auto tractIt = tract.begin(); std::vector profile; while(tractIt != tract.end()) { PointType p = *tractIt; // Get value from image profile.push_back( (mitk::ScalarType) imAccess.GetPixelByWorldCoordinates(p) ); ++tractIt; } profiles.push_back(profile); std::cout << std::endl; ++it; } if(profiles.size() == 0) return; m_IndividualProfiles = profiles; std::string title = "Fiber bundle plot"; this->SetPlotTitle( title.c_str() ); // initialize average profile std::vector averageProfile; std::vector profile = profiles.at(0); // can do this because we checked the size of profiles before for(unsigned int i=0; i profile = *profit; std::vector xAxis; for(unsigned int i=0; iInsertCurve( "" ); this->SetCurveData( curveId, xAxis, profile ); ++profit; id++; } m_Average = averageProfile; if(avg) { // Draw the average profile std::vector xAxis; for(unsigned int i=0; iInsertCurve( "" ); this->SetCurveData( curveId, xAxis, averageProfile ); QPen pen( Qt::SolidLine ); pen.setWidth(3); pen.setColor(Qt::red); this->SetCurvePen( curveId, pen ); id++; } this->Replot(); } 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/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTbssSkeletonizationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTbssSkeletonizationView.cpp index 3154134487..a15cd277a9 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTbssSkeletonizationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTbssSkeletonizationView.cpp @@ -1,502 +1,502 @@ /*=================================================================== 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. ===================================================================*/ // Qmitk #include "QmitkTbssSkeletonizationView.h" #include #include #include #include #include // mitk #include #include #include #include // Qt #include #include //vtk #include #include // Boost #include const std::string QmitkTbssSkeletonizationView::VIEW_ID = "org.mitk.views.tbssskeletonization"; using namespace berry; QmitkTbssSkeletonizationView::QmitkTbssSkeletonizationView() : QmitkAbstractView() , m_Controls( 0 ) , m_Activated(false) { } QmitkTbssSkeletonizationView::~QmitkTbssSkeletonizationView() { } void QmitkTbssSkeletonizationView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { //datamanager selection changed if (!this->IsActivated()) return; bool found3dImage = false; bool found4dImage = false; this->m_Controls->m_TubularName->setText(QString("Tubular Structure Mask: ")); this->m_Controls->m_TubularName->setEnabled(false); this->m_Controls->m_MeanLabel->setText(QString("Mean: ")); this->m_Controls->m_MeanLabel->setEnabled(false); this->m_Controls->m_PatientDataLabel->setText(QString("Patient Data: ")); this->m_Controls->m_PatientDataLabel->setEnabled(false); // iterate selection for (mitk::DataNode::Pointer node: nodes) { // only look at interesting types from valid nodes mitk::BaseData* nodeData = node->GetData(); std::string name = ""; node->GetStringProperty("name", name); if(nodeData) { if(QString("Image").compare(nodeData->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodeData); if(img->GetDimension() == 3) { bool isBinary(false); node->GetBoolProperty("binary", isBinary); if(isBinary) { QString label("Tubular Structure Mask: "); label.append(QString(name.c_str())); this->m_Controls->m_TubularName->setText(label); this->m_Controls->m_TubularName->setEnabled(true); } else { found3dImage = true; QString label("Mean: "); label.append(QString(name.c_str())); this->m_Controls->m_MeanLabel->setText(label); this->m_Controls->m_MeanLabel->setEnabled(true); } } else if(img->GetDimension() == 4) { found4dImage = true; QString label("Patient Data: "); label.append(QString(name.c_str())); this->m_Controls->m_PatientDataLabel->setText(label); this->m_Controls->m_PatientDataLabel->setEnabled(true); } } } } this->m_Controls->m_Skeletonize->setEnabled(found3dImage); this->m_Controls->m_Project->setEnabled(found3dImage && found4dImage); this->m_Controls->m_OutputMask->setEnabled(found3dImage && found4dImage); this->m_Controls->m_OutputDistanceMap->setEnabled(found3dImage && found4dImage); } void QmitkTbssSkeletonizationView::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::QmitkTbssSkeletonizationViewControls; m_Controls->setupUi( parent ); this->CreateConnections(); } } void QmitkTbssSkeletonizationView::SetFocus() { m_Controls->m_Skeletonize->setFocus(); } void QmitkTbssSkeletonizationView::Activated() { m_Activated = true; } void QmitkTbssSkeletonizationView::Deactivated() { m_Activated = false; } bool QmitkTbssSkeletonizationView::IsActivated() const { return m_Activated; } void QmitkTbssSkeletonizationView::Visible() { } void QmitkTbssSkeletonizationView::Hidden() { } void QmitkTbssSkeletonizationView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_Skeletonize), SIGNAL(clicked()), this, SLOT(Skeletonize() )); connect( (QObject*)(m_Controls->m_Project), SIGNAL(clicked()), this, SLOT(Project() )); } } void QmitkTbssSkeletonizationView::Skeletonize() { typedef itk::SkeletonizationFilter SkeletonisationFilterType; SkeletonisationFilterType::Pointer skeletonizer = SkeletonisationFilterType::New(); QList nodes = this->GetDataManagerSelection(); mitk::Image::Pointer meanImage = mitk::Image::New(); std::string name = ""; for (auto node: nodes) { // process only on valid nodes mitk::BaseData* nodeData = node->GetData(); if(nodeData) { if(QString("Image").compare(nodeData->GetNameOfClass())==0) { bool isBinary(false); node->GetBoolProperty("binary", isBinary); mitk::Image* img = static_cast(nodeData); if(img->GetDimension() == 3 && !isBinary) { meanImage = img; name = node->GetName(); } } } } // Calculate skeleton FloatImageType::Pointer itkImg = FloatImageType::New(); mitk::CastToItkImage(meanImage, itkImg); skeletonizer->SetInput(itkImg); skeletonizer->Update(); FloatImageType::Pointer output = skeletonizer->GetOutput(); mitk::Image::Pointer mitkOutput = mitk::Image::New(); mitk::CastToMitkImage(output, mitkOutput); name += "_skeleton"; AddToDataStorage(mitkOutput, name); } void QmitkTbssSkeletonizationView::Project() { typedef itk::SkeletonizationFilter SkeletonisationFilterType; typedef itk::ProjectionFilter ProjectionFilterType; typedef itk::DistanceMapFilter DistanceMapFilterType; SkeletonisationFilterType::Pointer skeletonizer = SkeletonisationFilterType::New(); QList nodes = this->GetDataManagerSelection(); mitk::Image::Pointer meanImage = mitk::Image::New(); mitk::Image::Pointer subjects = mitk::Image::New(); mitk::Image::Pointer tubular = mitk::Image::New(); for (auto node: nodes) { // process only on valid nodes mitk::BaseData* nodeData = node->GetData(); if(nodeData) { if(QString("Image").compare(nodeData->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodeData); if(img->GetDimension() == 3) { bool isBinary(false); node->GetBoolProperty("binary", isBinary); if(isBinary) { tubular = img; } else { meanImage = img; } } else if(img->GetDimension() == 4) { subjects = img; } } } } Float4DImageType::Pointer allFA; mitkPixelTypeMultiplex2(ConvertToItk,subjects->GetChannelDescriptor().GetPixelType(),subjects, allFA); // Calculate skeleton FloatImageType::Pointer itkImg = FloatImageType::New(); mitk::CastToItkImage(meanImage, itkImg); skeletonizer->SetInput(itkImg); skeletonizer->Update(); FloatImageType::Pointer output = skeletonizer->GetOutput(); mitk::Image::Pointer mitkOutput = mitk::Image::New(); mitk::CastToMitkImage(output, mitkOutput); AddToDataStorage(mitkOutput, "mean_FA_skeletonised"); // Retrieve direction image needed later by the projection filter DirectionImageType::Pointer directionImg = skeletonizer->GetVectorImage(); // Calculate distance image DistanceMapFilterType::Pointer distanceMapFilter = DistanceMapFilterType::New(); distanceMapFilter->SetInput(output); distanceMapFilter->Update(); FloatImageType::Pointer distanceMap = distanceMapFilter->GetOutput(); if(m_Controls->m_OutputDistanceMap->isChecked()) { mitk::Image::Pointer mitkDistance = mitk::Image::New(); mitk::CastToMitkImage(distanceMap, mitkDistance); AddToDataStorage(mitkDistance, "distance map"); } // Do projection // Ask a threshold to create a skeleton mask 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; ThresholdFilterType::Pointer thresholder = ThresholdFilterType::New(); thresholder->SetInput(output); thresholder->SetLowerThreshold(threshold); thresholder->SetUpperThreshold(std::numeric_limits::max()); thresholder->SetOutsideValue(0); thresholder->SetInsideValue(1); thresholder->Update(); CharImageType::Pointer thresholdedImg = thresholder->GetOutput(); if(m_Controls->m_OutputMask->isChecked()) { mitk::Image::Pointer mitkThresholded = mitk::Image::New(); mitk::CastToMitkImage(thresholdedImg, mitkThresholded); std::string maskName = "skeleton_mask_at_" + boost::lexical_cast(threshold); AddToDataStorage(mitkThresholded, maskName); } CharImageType::Pointer itkTubular = CharImageType::New(); mitk::CastToItkImage(tubular, itkTubular); ProjectionFilterType::Pointer projectionFilter = ProjectionFilterType::New(); projectionFilter->SetDistanceMap(distanceMap); projectionFilter->SetDirections(directionImg); projectionFilter->SetAllFA(allFA); projectionFilter->SetTube(itkTubular); projectionFilter->SetSkeleton(thresholdedImg); projectionFilter->Project(); Float4DImageType::Pointer projected = projectionFilter->GetProjections(); mitk::Image::Pointer mitkProjections = mitk::Image::New(); mitk::CastToMitkImage(projected, mitkProjections); AddToDataStorage(mitkProjections, "all_FA_projected"); } void QmitkTbssSkeletonizationView::AddToDataStorage(mitk::Image* img, std::string name) { mitk::DataNode::Pointer result = mitk::DataNode::New(); result->SetProperty( "name", mitk::StringProperty::New(name) ); result->SetData( img ); // add new image to data storage and set as active to ease further processing GetDataStorage()->Add( result ); } template -void QmitkTbssSkeletonizationView::ConvertToItk(mitk::PixelType ptype, mitk::Image* image, Float4DImageType::Pointer output) +void QmitkTbssSkeletonizationView::ConvertToItk(mitk::PixelType, mitk::Image* image, Float4DImageType::Pointer output) { output = Float4DImageType::New(); mitk::BaseGeometry* geo = image->GetGeometry(); mitk::Vector3D mitkSpacing = geo->GetSpacing(); mitk::Point3D mitkOrigin = geo->GetOrigin(); Float4DImageType::SpacingType spacing; spacing[0] = mitkSpacing[0]; spacing[1] = mitkSpacing[1]; spacing[2] = mitkSpacing[2]; spacing[3] = 1.0; // todo: check if spacing has length 4 Float4DImageType::PointType origin; origin[0] = mitkOrigin[0]; origin[1] = mitkOrigin[1]; origin[2] = mitkOrigin[2]; origin[3] = 0; Float4DImageType::SizeType size; size[0] = image->GetDimension(0); size[1] = image->GetDimension(1); size[2] = image->GetDimension(2); size[3] = image->GetDimension(3); Float4DImageType::DirectionType dir; vtkLinearTransform* lin = geo->GetVtkTransform(); vtkMatrix4x4 *m = lin->GetMatrix(); dir.Fill(0.0); for(int x=0; x<3; x++) { for(int y=0; y<3; y++) { dir[x][y] = m->GetElement(x,y); } } dir[3][3] = 1; output->SetSpacing(spacing); output->SetOrigin(origin); output->SetRegions(size); output->SetDirection(dir); output->Allocate(); if(image->GetDimension() == 4) { - int timesteps = image->GetDimension(3); + unsigned int timesteps = image->GetDimension(3); try{ // REPLACE THIS METHODE()ConvertToItk) WITH mitk::CastToItk // iterate through the subjects and copy data to output - for(int t=0; t inAcc(image,image->GetVolumeData(t)); - for(int x=0; xGetDimension(0); x++) + for(unsigned int x=0; xGetDimension(0); x++) { - for(int y=0; yGetDimension(1); y++) + for(unsigned int y=0; yGetDimension(1); y++) { - for(int z=0; zGetDimension(2); z++) + for(unsigned int z=0; zGetDimension(2); z++) { itk::Index<3> ix = {x, y, z}; itk::Index<4> ix4 = {x, y, z, t}; output->SetPixel(ix4, inAcc.GetPixelByIndex(ix)); } } } } } - catch(std::exception & e) + catch(const std::exception & e) { MITK_INFO << e.what(); } } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp index b2368f1c00..c263290d28 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging.tbss/src/internal/QmitkTractbasedSpatialStatisticsView.cpp @@ -1,1141 +1,1141 @@ /*=================================================================== 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. ===================================================================*/ // Qmitk #include "QmitkTractbasedSpatialStatisticsView.h" #include "mitkDataNodeObject.h" #include // Qt #include #include #include #include #include #include #include #include #include "vtkPoints.h" #include #include const std::string QmitkTractbasedSpatialStatisticsView::VIEW_ID = "org.mitk.views.tractbasedspatialstatistics"; using namespace berry; QmitkTractbasedSpatialStatisticsView::QmitkTractbasedSpatialStatisticsView() : QmitkAbstractView() , m_Controls( 0 ) , m_Activated(false) { } QmitkTractbasedSpatialStatisticsView::~QmitkTractbasedSpatialStatisticsView() { } void QmitkTractbasedSpatialStatisticsView::PerformChange() { m_Controls->m_RoiPlotWidget->ModifyPlot(m_Controls->m_Segments->value(), m_Controls->m_Average->isChecked()); } void QmitkTractbasedSpatialStatisticsView::OnSelectionChanged(berry::IWorkbenchPart::Pointer /*part*/, const QList& nodes) { //datamanager selection changed if (!this->IsActivated()) return; // Check which datatypes are selected in the datamanager and enable/disable widgets accordingly 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::FiberBundle* fib; - mitk::DataNode* start; - mitk::DataNode* end; + mitk::TbssRoiImage* roiImage = nullptr; + mitk::TbssImage* image = nullptr; + mitk::Image* img = nullptr; + mitk::FiberBundle* fib = nullptr; + mitk::DataNode* start = nullptr; + mitk::DataNode* end = nullptr; m_CurrentStartRoi = nullptr; m_CurrentEndRoi = nullptr; for (mitk::DataNode::Pointer node: nodes) { // only look at interesting types // check for valid data mitk::BaseData* nodeData = node->GetData(); 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("FiberBundle").compare(nodeData->GetNameOfClass())==0) { foundFiberBundle = true; fib = static_cast(nodeData); this->m_CurrentFiberNode = node; } if(QString("PlanarCircle").compare(nodeData->GetNameOfClass())==0) { if(!foundStartRoi) { start = node; this->m_CurrentStartRoi = node; foundStartRoi = true; } else { end = node; this->m_CurrentEndRoi = node; foundEndRoi = true; } } } } this->m_Controls->m_CreateRoi->setEnabled(found3dImage); this->m_Controls->m_ImportFsl->setEnabled(found4dImage); if(foundTbss && foundTbssRoi) { this->Plot(image, roiImage); } else if(found3dImage && foundFiberBundle && foundStartRoi && foundEndRoi) { this->PlotFiberBundle(fib, img, start, end); } else if(found3dImage && foundFiberBundle) { this->PlotFiberBundle(fib, img); } else if(foundTbss && foundStartRoi && foundEndRoi && foundFiberBundle) { this->PlotFiber4D(image, fib, start, end); } if(found3dImage) { this->InitPointsets(); } this->m_Controls->m_Cut->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi); this->m_Controls->m_SegmentLabel->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && (found3dImage || foundTbss)); this->m_Controls->m_Segments->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && (found3dImage || foundTbss)); this->m_Controls->m_Average->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); } void QmitkTractbasedSpatialStatisticsView::InitPointsets() { // Check if PointSetStart exsits, if not create it. m_P1 = this->GetDataStorage()->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->GetDataStorage()->Add(m_P1); m_Controls->m_PointWidget->SetPointSetNode(m_P1); auto renderWindowPart = this->GetRenderWindowPart(OPEN); auto axialSnc = renderWindowPart->GetQmitkRenderWindow("axial")->GetSliceNavigationController(); auto sagittalSnc = renderWindowPart->GetQmitkRenderWindow("sagittal")->GetSliceNavigationController(); auto coronalSnc = renderWindowPart->GetQmitkRenderWindow("coronal")->GetSliceNavigationController(); m_Controls->m_PointWidget->AddSliceNavigationController(axialSnc); m_Controls->m_PointWidget->AddSliceNavigationController(sagittalSnc); m_Controls->m_PointWidget->AddSliceNavigationController(coronalSnc); } } 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(); } // Table for the FSL TBSS import m_GroupModel = new QmitkTbssTableModel(); m_Controls->m_GroupInfo->setModel(m_GroupModel); } void QmitkTractbasedSpatialStatisticsView::SetFocus() { m_Controls->m_AddGroup->setFocus(); } void QmitkTractbasedSpatialStatisticsView::Activated() { m_Activated = true; } void QmitkTractbasedSpatialStatisticsView::Deactivated() { m_Activated = false; } bool QmitkTractbasedSpatialStatisticsView::IsActivated() const { return m_Activated; } void QmitkTractbasedSpatialStatisticsView::Visible() { } void QmitkTractbasedSpatialStatisticsView::Hidden() { } 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 QPointF&)), SLOT(Clicked(const QPointF&) ) ); connect( m_Controls->m_RoiPlotWidget->m_PlotPicker, SIGNAL(moved(const QPointF&)), SLOT(Clicked(const QPointF&) ) ); 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); QItemSelectionModel *selectionModel = temp->selectionModel(); QModelIndexList indices = selectionModel->selectedRows(); QModelIndex index; foreach(index, indices) { int row = index.row(); m_GroupModel->removeRows(row, 1, QModelIndex()); } } 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); } } 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 = ""; QList nodes = this->GetDataManagerSelection(); for (auto node: nodes) { 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 GetDataStorage()->Add( result ); // show the results mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } void QmitkTractbasedSpatialStatisticsView::Clicked(const QPointF& pos) { int index = (int)pos.x(); if(m_Roi.size() > 0 && m_CurrentGeometry != nullptr && !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]; this->GetRenderWindowPart()->SetSelectedPosition(p); m_Controls->m_RoiPlotWidget->drawBar(index); } else if(m_Controls->m_RoiPlotWidget->IsPlottingFiber() ) { mitk::Point3D point = m_Controls->m_RoiPlotWidget->GetPositionInWorld(index); this->GetRenderWindowPart()->SetSelectedPosition(point); } } void QmitkTractbasedSpatialStatisticsView::Cut() { mitk::BaseData* fibData = m_CurrentFiberNode->GetData(); mitk::FiberBundle* fib = static_cast(fibData); mitk::PlaneGeometry* startGeometry2D = const_cast(dynamic_cast(m_CurrentStartRoi->GetData())->GetPlaneGeometry()); mitk::PlaneGeometry* endGeometry2D = const_cast(dynamic_cast(m_CurrentEndRoi->GetData())->GetPlaneGeometry()); mitk::Point3D startCenter = dynamic_cast(m_CurrentStartRoi->GetData())->GetWorldControlPoint(0); //center Point of start roi mitk::Point3D endCenter = dynamic_cast(m_CurrentEndRoi->GetData())->GetWorldControlPoint(0); //center Point of end roi mitk::FiberBundle::Pointer inStart = fib->ExtractFiberSubset(m_CurrentStartRoi, nullptr); mitk::FiberBundle::Pointer inBoth = inStart->ExtractFiberSubset(m_CurrentEndRoi, nullptr); 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(nullptr); 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::FiberBundle::Pointer cutBundle = mitk::FiberBundle::New(polyData); mitk::DataNode::Pointer cutNode = mitk::DataNode::New(); cutNode->SetData(cutBundle); std::string name = "fiberCut"; cutNode->SetName(name); GetDataStorage()->Add(cutNode); } void QmitkTractbasedSpatialStatisticsView::CreateRoi() { 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; QList nodes = this->GetDataManagerSelection(); for (auto node: nodes) { 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); m_PointSetNode = this->m_Controls->m_PointWidget->GetPointSet(); // Set Pointset to analyzer analyzer.SetPointSet(m_PointSetNode); // Run Analyzer try { analyzer.MakeRoi(); } catch (const mitk::Exception& e) { QMessageBox msgBox; msgBox.setText(QString::fromStdString(e.what())); msgBox.exec(); } // 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::PlotFiber4D(mitk::TbssImage* image, mitk::FiberBundle* fib, mitk::DataNode* startRoi, mitk::DataNode* endRoi) { if(m_Controls->m_TabWidget->currentWidget() == m_Controls->m_MeasureTAB) { m_CurrentGeometry = image->GetGeometry(); m_Controls->m_RoiPlotWidget->SetGroups(image->GetGroupInfo()); m_Controls->m_RoiPlotWidget->SetProjections(image->GetImage()); m_Controls->m_RoiPlotWidget->SetMeasure( image->GetMeasurementInfo() ); m_Controls->m_RoiPlotWidget->PlotFiber4D(image, fib, startRoi, endRoi, m_Controls->m_Segments->value()); } } void QmitkTractbasedSpatialStatisticsView:: PlotFiberBundle(mitk::FiberBundle *fib, mitk::Image* img, mitk::DataNode* startRoi, mitk::DataNode* 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 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(); } m_Controls->m_RoiPlotWidget->SetPlottingFiber(false); }