diff --git a/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp b/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp index c7ab1ff07e..a790432044 100644 --- a/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp +++ b/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp @@ -1,119 +1,160 @@ /*=================================================================== 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 __mitkNrrdTbssImageWriter__cpp #define __mitkNrrdTbssImageWriter__cpp #include "mitkNrrdTbssImageWriter.h" #include "itkMetaDataDictionary.h" #include "itkMetaDataObject.h" #include "itkNrrdImageIO.h" //#include "itkNiftiImageIO.h" #include "itkImageFileWriter.h" #include "itksys/SystemTools.hxx" #include "boost/lexical_cast.hpp" #include #include mitk::NrrdTbssImageWriter::NrrdTbssImageWriter() : m_FileName(""), m_FilePrefix(""), m_FilePattern(""), m_Success(false) { this->SetNumberOfRequiredInputs( 1 ); } mitk::NrrdTbssImageWriter::~NrrdTbssImageWriter() {} void mitk::NrrdTbssImageWriter::GenerateData() { m_Success = false; InputType* input = this->GetInput(); if (input == NULL) { itkWarningMacro(<<"Sorry, input to NrrdTbssImageWriter is NULL!"); return; } if ( m_FileName == "" ) { itkWarningMacro( << "Sorry, filename has not been set!" ); return ; } itk::VectorImage::Pointer img = input->GetImage(); + + + std::string key; + std::string val; + + + /* For the case of a tbss image containing data of the patients: + Save info about the groups and the type of measurement */ + + + + + std::vector< std::pair > groups = input->GetGroupInfo(); + std::vector< std::pair >::iterator it = groups.begin(); + + int i=0; + while(it != groups.end()) + { + std::pair p = *it; + + key = "Group_index_" + boost::lexical_cast(i); + val = " " + p.first + " " + boost::lexical_cast(p.second); + //sprintf( keybuffer, "Group_index_%04d", std::string(i) ); + // sprintf( valbuffer, "%1d %1d", p.first, p.second); + + //std::cout << valbuffer << std::endl; + + //itk::EncapsulateMetaData< std::string >(input->GetImage()->GetMetaDataDictionary(),std::string(keybuffer),std::string(valbuffer)); + itk::EncapsulateMetaData< std::string >(input->GetImage()->GetMetaDataDictionary(),key,val); + it++; + ++i; + } + + key = "Measurement info"; + val = input->GetMeasurementInfo(); + itk::EncapsulateMetaData< std::string >(input->GetImage()->GetMetaDataDictionary(),key,val); + + + + + typedef itk::VectorImage ImageType; itk::NrrdImageIO::Pointer io = itk::NrrdImageIO::New(); io->SetFileType( itk::ImageIOBase::Binary ); io->UseCompressionOn(); typedef itk::ImageFileWriter WriterType; WriterType::Pointer nrrdWriter = WriterType::New(); nrrdWriter->UseInputMetaDataDictionaryOn(); nrrdWriter->SetInput( img ); nrrdWriter->SetImageIO(io); nrrdWriter->SetFileName(m_FileName); nrrdWriter->UseCompressionOn(); nrrdWriter->SetImageIO(io); try { nrrdWriter->Update(); } catch (itk::ExceptionObject e) { std::cout << e << std::endl; } m_Success = true; } void mitk::NrrdTbssImageWriter::SetInput( InputType* tbssVol ) { this->ProcessObject::SetNthInput( 0, tbssVol ); } mitk::TbssImage* mitk::NrrdTbssImageWriter::GetInput() { if ( this->GetNumberOfInputs() < 1 ) { return NULL; } else { return dynamic_cast ( this->ProcessObject::GetInput( 0 ) ); } } std::vector mitk::NrrdTbssImageWriter::GetPossibleFileExtensions() { std::vector possibleFileExtensions; possibleFileExtensions.push_back(".tbss"); return possibleFileExtensions; } #endif //__mitkNrrdTbssImageWriter__cpp diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp index 118e7e599c..f56b5e8eac 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp @@ -1,968 +1,980 @@ /*=================================================================== 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 #include #include #include #include #include #include #include #include #include #include #include + + + QmitkTbssRoiAnalysisWidget::QmitkTbssRoiAnalysisWidget( QWidget * parent ) : QmitkPlotWidget(parent) { m_PlotPicker = new QwtPlotPicker(m_Plot->canvas()); m_PlotPicker->setSelectionFlags(QwtPicker::PointSelection | QwtPicker::ClickSelection | QwtPicker::DragSelection); m_PlotPicker->setTrackerMode(QwtPicker::ActiveOnly); m_PlottingFiberBundle = false; } void QmitkTbssRoiAnalysisWidget::DoPlotFiberBundles(mitk::FiberBundleX *fib, mitk::Image* img, mitk::PlanarFigure* startRoi, mitk::PlanarFigure* endRoi, bool avg, int number) { - mitk::Geometry3D* currentGeometry = fib->GetGeometry(); + TractContainerType tracts = CreateTracts(fib, startRoi, endRoi); - mitk::PlaneGeometry* startGeometry2D = dynamic_cast( const_cast(startRoi->GetGeometry2D()) ); - mitk::PlaneGeometry* endGeometry2D = dynamic_cast( const_cast(endRoi->GetGeometry2D()) ); + TractContainerType resampledTracts = ParameterizeTracts(tracts, number); + // Now we have the resampled tracts. Next we should use these points to read out the values - 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); + PlotFiberBundles(resampledTracts, img, avg); + m_CurrentTracts = resampledTracts; +} - int num = inBoth->GetNumFibers(); +TractContainerType QmitkTbssRoiAnalysisWidget::CreateTracts(mitk::FiberBundleX *fib, + mitk::PlanarFigure *startRoi, + mitk::PlanarFigure *endRoi) +{ + mitk::PlaneGeometry* startGeometry2D = dynamic_cast( const_cast(startRoi->GetGeometry2D()) ); + mitk::PlaneGeometry* endGeometry2D = dynamic_cast( const_cast(endRoi->GetGeometry2D()) ); - TractContainerType tracts; + 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); - vtkSmartPointer fiberPolyData = inBoth->GetFiberPolyData(); - vtkCellArray* lines = fiberPolyData->GetLines(); - lines->InitTraversal(); + int num = inBoth->GetNumFibers(); - // 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(NULL); - lines->GetNextCell ( numPointsInCell, pointsInCell ); - int startId = 0; - int endId = 0; + TractContainerType tracts; - float minDistStart = std::numeric_limits::max(); - float minDistEnd = std::numeric_limits::max(); + vtkSmartPointer fiberPolyData = inBoth->GetFiberPolyData(); + vtkCellArray* lines = fiberPolyData->GetLines(); + lines->InitTraversal(); - for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) + + // 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(NULL); + lines->GetNextCell ( numPointsInCell, pointsInCell ); + int startId = 0; + int endId = 0; - double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); + float minDistStart = std::numeric_limits::max(); + float minDistEnd = std::numeric_limits::max(); - 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) + for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { - minDistStart = distanceToStart; - startId = pointInCellID; - } - if(distanceToEnd < minDistEnd) - { - minDistEnd = distanceToEnd; - endId = 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); - /* 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(distanceToStart < minDistStart) + { + minDistStart = distanceToStart; + startId = pointInCellID; + } - TractType singleTract; - PointType point; + if(distanceToEnd < minDistEnd) + { + minDistEnd = distanceToEnd; + endId = pointInCellID; + } - 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 - 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(); + /* 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*/ - mitk::Point3D pStart; - pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; + TractType singleTract; + PointType point; - mitk::Point3D pSecond; - pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; + if(startId < endId) + { - bool startOnPositive = startGeometry2D->IsAbove(pStart); - bool secondOnPositive = startGeometry2D->IsAbove(pSecond); + // Calculate the intersection of the ROI with the startRoi and decide if the startId is part of the roi or must be cut of + 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]; - mitk::Vector3D onPlane; - onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[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]; - 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]; - } + bool startOnPositive = startGeometry2D->IsAbove(pStart); + bool secondOnPositive = startGeometry2D->IsAbove(pSecond); - double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); - mitk::Vector3D newPoint = (p0-p1); + mitk::Vector3D onPlane; + onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; - 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*/ + 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 *start = fiberPolyData->GetPoint( pointsInCell[startId] ); - point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; - singleTract.push_back(point); - } + double 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]; - for( int pointInCellID( startId+1 ); pointInCellID < endId ; pointInCellID++) - { - // push back point - double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); - point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; - singleTract.push_back( point ); + 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*/ - /* endId must be included if endId and endId-1 lie on the same side of the - plane defined by endRoi*/ + double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); + point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; + singleTract.push_back(point); + } - 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]; + for( int pointInCellID( startId+1 ); pointInCellID < endId ; pointInCellID++) + { + // push back point + double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); + point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; + singleTract.push_back( point ); - 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]; + /* endId must be included if endId and endId-1 lie on the same side of the + plane defined by endRoi*/ - normal = endGeometry2D->GetNormal(); - bool lastOnPositive = endGeometry2D->IsAbove(pLast); - bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); + p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); + p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; - onPlane[0] = endCenter[0]; - onPlane[1] = endCenter[1]; - onPlane[2] = endCenter[2]; + p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); + p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[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 */ + mitk::Point3D pLast; + pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; - p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); - point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; - singleTract.push_back( point ); + mitk::Point3D pBeforeLast; + pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; - p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); - } + normal = endGeometry2D->GetNormal(); + bool lastOnPositive = endGeometry2D->IsAbove(pLast); + bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); - d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); - newPoint = (p0-p1); + onPlane[0] = endCenter[0]; + onPlane[1] = endCenter[1]; + onPlane[2] = endCenter[2]; - point[0] = d*newPoint[0] + p0[0]; - point[1] = d*newPoint[1] + p0[1]; - point[2] = d*newPoint[2] + p0[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 */ - singleTract.push_back(point); + 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 ); - } - else{ + newPoint = (p0-p1); - // Calculate the intersection of the ROI with the startRoi and decide if the startId is part of the roi or must be cut of - double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); - mitk::Vector3D p0; - p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; + point[0] = d*newPoint[0] + p0[0]; + point[1] = d*newPoint[1] + p0[1]; + point[2] = d*newPoint[2] + p0[2]; - p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); - mitk::Vector3D p1; - p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; + singleTract.push_back(point); - // 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]; + } + else{ - bool startOnPositive = startGeometry2D->IsAbove(pStart); - bool secondOnPositive = startGeometry2D->IsAbove(pSecond); + // Calculate the intersection of the ROI with the startRoi and decide if the startId is part of the roi or must be cut of + double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); + mitk::Vector3D p0; + p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; - mitk::Vector3D onPlane; - onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[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(); - 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::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]; - double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); - mitk::Vector3D newPoint = (p0-p1); + 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]; - point[0] = d*newPoint[0] + p0[0]; - point[1] = d*newPoint[1] + p0[1]; - point[2] = d*newPoint[2] + p0[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]; + } - 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*/ + double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); + mitk::Vector3D newPoint = (p0-p1); - double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); - point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; - singleTract.push_back(point); - } + point[0] = d*newPoint[0] + p0[0]; + point[1] = d*newPoint[1] + p0[1]; + point[2] = d*newPoint[2] + p0[2]; - for( int pointInCellID( startId-1 ); pointInCellID > endId ; pointInCellID--) - { - // push back point - double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); - point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; - singleTract.push_back( point ); + 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*/ - /* endId must be included if endId and endI+1 lie on the same side of the - plane defined by endRoi*/ + double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); + point[0] = start[0]; point[1] = start[1]; point[2] = start[2]; + singleTract.push_back(point); + } - 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]; + for( int pointInCellID( startId-1 ); pointInCellID > endId ; pointInCellID--) + { + // push back point + double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); + point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; + singleTract.push_back( point ); - 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]; + /* endId must be included if endId and endI+1 lie on the same side of the + plane defined by endRoi*/ - bool lastOnPositive = endGeometry2D->IsAbove(pLast); - bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); - normal = endGeometry2D->GetNormal(); + 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]; - onPlane[0] = endCenter[0]; - onPlane[1] = endCenter[1]; - onPlane[2] = endCenter[2]; + mitk::Point3D pLast; + pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[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 */ + mitk::Point3D pBeforeLast; + pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; - p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); - point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; - singleTract.push_back( point ); + bool lastOnPositive = endGeometry2D->IsAbove(pLast); + bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); + normal = endGeometry2D->GetNormal(); - p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); - } - d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); - newPoint = (p0-p1); + onPlane[0] = endCenter[0]; + onPlane[1] = endCenter[1]; + onPlane[2] = endCenter[2]; - point[0] = d*newPoint[0] + p0[0]; - point[1] = d*newPoint[1] + p0[1]; - point[2] = d*newPoint[2] + p0[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 */ - singleTract.push_back(point); + 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 ); - tracts.push_back(singleTract); + 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); + } - //todo: Make number of samples selectable by user - TractContainerType resampledTracts = ParameterizeTracts(tracts, number); - // Now we have the resampled tracts. Next we should use these points to read out the values + tracts.push_back(singleTract); + } - PlotFiberBundles(resampledTracts, img, avg); - m_CurrentTracts = resampledTracts; + return tracts; } void QmitkTbssRoiAnalysisWidget::PlotFiberBetweenRois(mitk::FiberBundleX *fib, mitk::Image* img, mitk::PlanarFigure* startRoi, mitk::PlanarFigure* endRoi, bool avg, int number) { if(fib == NULL || img == NULL || startRoi == NULL || endRoi == NULL) 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 == NULL || m_CurrentImage == NULL || m_CurrentStartRoi == NULL || m_CurrentEndRoi == NULL) + if(m_Fib == NULL || m_CurrentTbssImage == NULL || m_CurrentStartRoi == NULL || m_CurrentEndRoi == NULL) return; - DoPlotFiberBundles(m_Fib, m_CurrentImage, m_CurrentStartRoi, m_CurrentEndRoi, avg, number); - - + 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(TractContainerType::iterator it = tracts.begin(); it != tracts.end(); ++it) { TractType resampledTract; TractType tract = *it; // Calculate the total length float totalLength = 0; if(tract.size() < 2) continue; PointType p0 = tract.at(0); for(int i = 1; i 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) { - TractContainerType tractsAtIndex; + float xSum = 0.0; float ySum = 0.0; float zSum = 0.0; for(TractContainerType::iterator 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(); float xPos = xSum / number; float yPos = ySum / number; float zPos = zSum / number; mitk::Point3D pos; pos[0] = xPos; pos[1] = yPos; pos[2] = zPos; return pos; } -std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfiles(std::string preprocessed) +std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfiles() { MITK_INFO << "make profiles!"; std::vector< std::vector > profiles; - //No results were preprocessed, so they must be calculated now. - if(preprocessed == "") + + int size = m_Projections->GetVectorLength(); + for(int s=0; sGetVectorLength(); - for(int s=0; s profile; + RoiType::iterator it; + it = m_Roi.begin(); + while(it != m_Roi.end()) { - // Iterate trough the roi - std::vector profile; - RoiType::iterator it; - it = m_Roi.begin(); - while(it != m_Roi.end()) - { - itk::Index<3> ix = *it; - profile.push_back(m_Projections->GetPixel(ix).GetElement(s)); - it++; - } - int pSize = profile.size(); - profiles.push_back(profile); - } - } - else{ // Use preprocessed results - std::ifstream file(preprocessed.c_str()); - if(file.is_open()) - { - std::string line; - while(getline(file,line)) - { - std::vector tokens; - Tokenize(line, tokens); - std::vector::iterator it; - it = tokens.begin(); - std::vector< double > profile; - while(it != tokens.end()) - { - std::string s = *it; - profile.push_back (atof( s.c_str() ) ); - ++it; - } - profiles.push_back(profile); - } + 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 - int nprof = profiles.size(); - while(it != m_Groups.end() && profiles.size() > 0) { std::pair p = *it; int size = p.second; //initialize a vector of the right length with zeroes std::vector averageProfile; for(int i=0; iClear(); + std::vector > groupProfiles = CalculateGroupProfiles(); + Plot(groupProfiles); +} +void QmitkTbssRoiAnalysisWidget::Plot(std::vector > groupProfiles) +{ + this->Clear(); + m_Vals.clear(); + std::vector v1; - m_Vals.clear(); + std::vector xAxis; + for(int i=0; i v1; + // fill m_Vals. This might be used by the user to copy data to the clipboard + for(int i=0; i > groupProfiles = CalculateGroupProfiles(preprocessed); + v1.clear(); - std::vector xAxis; - for(int i=0; iSetPlotTitle( title.c_str() ); + QPen pen( Qt::SolidLine ); + pen.setWidth(2); - for(int j=0; j >::iterator it; + it = m_Groups.begin(); - } + int c = 0; //the current profile number - std::string title = m_Measure + " profiles on the "; - title.append(m_Structure); - this->SetPlotTitle( title.c_str() ); - QPen pen( Qt::SolidLine ); - pen.setWidth(2); + 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; - std::vector< std::pair >::iterator it; - it = m_Groups.begin(); + pen.setColor(colors[c]); + int curveId = this->InsertCurve( group.first.c_str() ); + this->SetCurveData( curveId, xAxis, groupProfiles.at(c) ); - int c = 0; //the current profile number - QColor colors[4] = {Qt::green, Qt::blue, Qt::yellow, Qt::red}; + this->SetCurvePen( curveId, pen ); - while(it != m_Groups.end() && groupProfiles.size() > 0) - { + c++; + it++; - 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) ); + QwtLegend *legend = new QwtLegend; + this->SetLegend(legend, QwtPlot::RightLegend, 0.5); - this->SetCurvePen( curveId, pen ); - c++; - it++; + std::cout << m_Measure << std::endl; + this->m_Plot->setAxisTitle(0, m_Measure.c_str()); + this->m_Plot->setAxisTitle(3, "Position"); - } + this->Replot(); +} - QwtLegend *legend = new QwtLegend; - this->SetLegend(legend, QwtPlot::RightLegend, 0.5); +std::vector< std::vector > QmitkTbssRoiAnalysisWidget::CalculateGroupProfilesFibers(mitk::TbssImage::Pointer tbssImage, + mitk::FiberBundleX *fib, + mitk::PlanarFigure* startRoi, + mitk::PlanarFigure* endRoi, + int number) +{ + TractContainerType tracts = CreateTracts(fib, startRoi, endRoi); + TractContainerType resampledTracts = ParameterizeTracts(tracts, number); - std::cout << m_Measure << std::endl; - this->m_Plot->setAxisTitle(0, m_Measure.c_str()); - this->m_Plot->setAxisTitle(3, "Position"); + int nTracts = resampledTracts.size(); - this->Replot(); + this->Clear(); -} + // For every group we have m fibers * n subjects of profiles to fill + std::vector< std::vector > profiles; -void QmitkTbssRoiAnalysisWidget::PlotFiberBundles(TractContainerType tracts, mitk::Image *img, bool avg) -{ - this->Clear(); - std::vector::iterator it = tracts.begin(); + // calculate individual profiles by going through all n subjects + int size = m_Projections->GetVectorLength(); + for(int s=0; s profile; + TractType::iterator 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)); - int min = std::numeric_limits::max(); - TractType smallestTract; - while(it != tracts.end()) - { - TractType tract = *it; - if(tract.size() correspondingIndices; - TractType correspondingPoints; + 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 - for(int i=0; i > groupProfiles; + + while(it != m_Groups.end() && profiles.size() > 0) { - PointType p = smallestTract[i]; + std::pair p = *it; + int size = p.second; + + //initialize a vector of the right length with zeroes + std::vector averageProfile; + for(int i=0; i::max(); - int correspondingIndex = 0; - PointType correspondingPoint; + // Average the right number of profiles - // Search for the point on the second tract with the smallest distance - // to p and memorize it - for(int j=0; j > groupProfiles = CalculateGroupProfilesFibers(tbssImage, fib, startRoi, endRoi, number); + Plot(groupProfiles); - for(int i=0; iClear(); + std::vector::iterator it = tracts.begin(); - // Get the values along the curves from a 3D images. should also contain info about the position on screen. std::vector< std::vector > profiles; it = tracts.begin(); while(it != tracts.end()) { - std::cout << "Tract\n"; + TractType tract = *it; TractType::iterator tractIt = tract.begin(); std::vector profile; while(tractIt != tract.end()) { PointType p = *tractIt; - std::cout << p[0] << ' ' << p[1] << ' ' << p[2] << '\n'; - // Get value from image profile.push_back( (double)img->GetPixelValueByWorldCoordinate(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(int i=0; i >::iterator profit = profiles.begin(); int id=0; while(profit != profiles.end()) { std::vector profile = *profit; std::vector xAxis; for(int i=0; iInsertCurve( QString::number(id).toStdString().c_str() ); + int curveId = this->InsertCurve( "" ); this->SetCurveData( curveId, xAxis, profile ); ++profit; id++; } m_Average = averageProfile; if(avg) { // Draw the average profile std::vector xAxis; for(int i=0; iInsertCurve( QString::number(id).toStdString().c_str() ); + int curveId = this->InsertCurve( "" ); 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/src/QmitkTbssRoiAnalysisWidget.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h index 6c85b4f9d7..cf283339cd 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h @@ -1,217 +1,231 @@ /*=================================================================== 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 QmitkTbssRoiAnalysisWidget_H_ #define QmitkTbssRoiAnalysisWidget_H_ #include "QmitkPlotWidget.h" #include #include "QmitkExtExports.h" #include "mitkPlanarFigure.h" #include "itkVectorImage.h" #include #include +#include #include #include typedef itk::VectorImage VectorImageType; typedef std::vector< itk::Index<3> > RoiType; typedef itk::Point PointType; typedef std::vector< PointType> TractType; typedef std::vector< TractType > TractContainerType; /** * \brief Plot widget for TBSS Data * This widget can plot regions of interest on TBSS projection data. */ class DIFFUSIONIMAGING_EXPORT QmitkTbssRoiAnalysisWidget : public QmitkPlotWidget { Q_OBJECT public: QmitkTbssRoiAnalysisWidget( QWidget * parent); virtual ~QmitkTbssRoiAnalysisWidget(); /* \brief Set group information */ void SetGroups(std::vector< std::pair > groups) { m_Groups = groups; } /* \brief Draws the group averaged profiles */ - void DrawProfiles(std::string preprocessed); + void DrawProfiles(); + + + void PlotFiber4D(mitk::TbssImage::Pointer tbssImage, + mitk::FiberBundleX *fib, + mitk::PlanarFigure* startRoi, + mitk::PlanarFigure* endRoi, + int number); void PlotFiberBundles(TractContainerType tracts, mitk::Image* img, bool avg=false); /* \brief Sets the projections of the individual subjects */ void SetProjections(VectorImageType::Pointer projections) { m_Projections = projections; } /* \brief Set the region of interest*/ void SetRoi(RoiType roi) { m_Roi = roi; } /* \brief Set structure information to display in the plot */ void SetStructure(std::string structure) { m_Structure = structure; } /* \brief Set measurement type for display in the plot */ void SetMeasure(std::string measure) { m_Measure = measure; } /* \brief Draws a bar to indicate were the user clicked in the plot */ void drawBar(int x); /* \brief Returns the values of the group averaged profiles */ std::vector > GetVals() { return m_Vals; } /* \brief Returns the values of the individual subjects profiles */ std::vector > GetIndividualProfiles() { return m_IndividualProfiles; } std::vector GetAverageProfile() { return m_Average; } void SetPlottingFiber(bool b) { m_PlottingFiberBundle = b; } bool IsPlottingFiber() { return m_PlottingFiberBundle; } void PlotFiberBetweenRois(mitk::FiberBundleX *fib, mitk::Image* img, mitk::PlanarFigure* startRoi, mitk::PlanarFigure* endRoi, bool avg=-1, int number=25); // Takes an index which is an x coordinate from the plot and finds the corresponding position in world space mitk::Point3D GetPositionInWorld(int index); void ModifyPlot(int number, bool avg); QwtPlotPicker* m_PlotPicker; protected: mitk::FiberBundleX* m_Fib; std::vector< std::vector > m_Vals; std::vector< std::vector > m_IndividualProfiles; std::vector< double > m_Average; - std::vector< std::vector > CalculateGroupProfiles(std::string preprocessed); + std::vector< std::vector > CalculateGroupProfiles(); + std::vector< std::vector > CalculateGroupProfilesFibers(mitk::TbssImage::Pointer tbssImage, + mitk::FiberBundleX *fib, + mitk::PlanarFigure* startRoi, + mitk::PlanarFigure* endRoi, + int number); + void Plot(std::vector > groupProfiles); void Tokenize(const std::string& str, std::vector& tokens, const std::string& delimiters = " ") { // Skip delimiters at beginning. std::string::size_type lastPos = str.find_first_not_of(delimiters, 0); // Find first "non-delimiter". std::string::size_type pos = str.find_first_of(delimiters, lastPos); while (std::string::npos != pos || std::string::npos != lastPos) { // Found a token, add it to the vector. tokens.push_back(str.substr(lastPos, pos - lastPos)); // Skip delimiters. Note the "not_of" lastPos = str.find_first_not_of(delimiters, pos); // Find next "non-delimiter" pos = str.find_first_of(delimiters, lastPos); } } std::vector< std::pair > m_Groups; VectorImageType::Pointer m_Projections; RoiType m_Roi; std::string m_Structure; std::string m_Measure; bool m_PlottingFiberBundle; // true when the plot results from a fiber tracking result (vtk .fib file) // Resample a collection of tracts so that every tract contains #number equidistant samples TractContainerType ParameterizeTracts(TractContainerType tracts, int number); - - TractContainerType m_CurrentTracts; - mitk::Image* m_CurrentImage; + mitk::TbssImage* m_CurrentTbssImage; mitk::PlanarFigure* m_CurrentStartRoi; mitk::PlanarFigure* m_CurrentEndRoi; void DoPlotFiberBundles(mitk::FiberBundleX *fib, mitk::Image* img, mitk::PlanarFigure* startRoi, mitk::PlanarFigure* endRoi, bool avg=false, int number=25); + /* \brief Creates tracts from a mitk::FiberBundleX and two planar figures indicating the start end end point */ + TractContainerType CreateTracts(mitk::FiberBundleX *fib, mitk::PlanarFigure* startRoi, mitk::PlanarFigure* endRoi); }; #endif 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 db6de5555f..389d500d13 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp @@ -1,1132 +1,1125 @@ /*=================================================================== 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 "QmitkStdMultiWidget.h" #include "mitkDataNodeObject.h" #include // Qt #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() : 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"); - } -*/ // 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::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) + else if(found3dImage && foundFiberBundle && foundStartRoi && foundEndRoi) { this->PlotFiberBundle(fib, img, start, end); } - else if(found3dImage == true && foundFiberBundle) + 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); - this->m_Controls->m_Segments->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); + 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->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(); } // Table for the FSL TBSS import m_GroupModel = new QmitkTbssTableModel(); m_Controls->m_GroupInfo->setModel(m_GroupModel); } void QmitkTractbasedSpatialStatisticsView::Activated() { QmitkFunctionality::Activated(); } 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); 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 = ""; 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::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; 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.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::PlotFiber4D(mitk::TbssImage* image, + mitk::FiberBundleX* fib, + mitk::PlanarFigure* startRoi, + mitk::PlanarFigure* 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::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->DrawProfiles(); } m_Controls->m_RoiPlotWidget->SetPlottingFiber(false); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h index 0f3a45cfff..c0c3e50dfd 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h @@ -1,170 +1,175 @@ /*=================================================================== 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 QmitkTractbasedSpatialStatisticsView_h #define QmitkTractbasedSpatialStatisticsView_h #include #include "ui_QmitkTractbasedSpatialStatisticsViewControls.h" #include #include #include #include #include #include #include #include "QmitkTbssTableModel.h" #include "QmitkTbssMetaTableModel.h" #include // Image types typedef short DiffusionPixelType; typedef itk::Image CharImageType; typedef itk::Image UCharImageType; typedef itk::Image Float4DImageType; typedef itk::Image FloatImageType; typedef itk::VectorImage VectorImageType; // Readers/Writers typedef itk::ImageFileReader< CharImageType > CharReaderType; typedef itk::ImageFileReader< UCharImageType > UCharReaderType; typedef itk::ImageFileWriter< CharImageType > CharWriterType; typedef itk::ImageFileReader< FloatImageType > FloatReaderType; typedef itk::ImageFileWriter< FloatImageType > FloatWriterType; typedef itk::ImageFileReader< Float4DImageType > Float4DReaderType; typedef itk::ImageFileWriter< Float4DImageType > Float4DWriterType; /*! * \brief This plugin provides an extension for Tract-based spatial statistics (see Smith et al., 2009. http://dx.doi.org/10.1016/j.neuroimage.2006.02.024) * TBSS enables analyzing the brain by a pipeline of registration, skeletonization, and projection that results in a white matter skeleton * for all subjects that are analyzed statistically. This plugin provides functionality to select single tracts and analyze them seperately. */ class QmitkTractbasedSpatialStatisticsView : public QmitkFunctionality { Q_OBJECT public: static const std::string VIEW_ID; QmitkTractbasedSpatialStatisticsView(); virtual ~QmitkTractbasedSpatialStatisticsView(); virtual void CreateQtPartControl(QWidget *parent); /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); /// \brief Called when the functionality is activated virtual void Activated(); virtual void Deactivated(); protected slots: // Creates Roi void CreateRoi(); void Clicked(const QwtDoublePoint& pos); // Import of FSL TBSS data void TbssImport(); // Add a group as metadata. This metadata is required by the plotting functionality void AddGroup(); // Remove a group void RemoveGroup(); // Copies the values displayed in the plot widget to clipboard, i.e. exports the data void CopyToClipboard(); // Method to cut away parts of fiber bundles that should not be plotted. void Cut(); // Adjust plot widget void PerformChange(); protected: /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); // Creates a plot using a 4D image containing the projections of all subjects and a region of interest void Plot(mitk::TbssImage*, mitk::TbssRoiImage*); void PlotFiberBundle(mitk::FiberBundleX* fib, mitk::Image* img, mitk::PlanarFigure* startRoi=NULL, mitk::PlanarFigure* endRoi=NULL); + void PlotFiber4D(mitk::TbssImage*, mitk::FiberBundleX* fib, mitk::PlanarFigure* startRoi=NULL, mitk::PlanarFigure* endRoi=NULL); + // Create a point set. This point set defines the points through which a region of interest should go void InitPointsets(); // Pointset and DataNode to contain the PointSet used in ROI creation mitk::PointSet::Pointer m_PointSetNode; mitk::DataNode::Pointer m_P1; // GUI widgets Ui::QmitkTractbasedSpatialStatisticsViewControls* m_Controls; /* A pointer to the QmitkStdMultiWidget. Used for interaction with the plot widget (clicking in the plot widget makes the image cross jump to the corresponding location on the skeleton).*/ QmitkStdMultiWidget* m_MultiWidget; // Used to save the region of interest in a vector of itk::index. std::vector< itk::Index<3> > m_Roi; mitk::FiberBundleX* m_Fib; mitk::Geometry3D* m_CurrentGeometry; // A table model for saving group information in a name,number pair. QmitkTbssTableModel* m_GroupModel; // Convenience function for adding a new image to the datastorage and giving it a name. void AddTbssToDataStorage(mitk::Image* image, std::string name); mitk::DataNode::Pointer m_CurrentFiberNode; // needed for the index property when interacting with the plot widget // needed when a plot should only show values between a start end end roi mitk::DataNode::Pointer m_CurrentStartRoi; mitk::DataNode::Pointer m_CurrentEndRoi; + + + }; #endif // _QMITKTRACTBASEDSPATIALSTATISTICSVIEW_H_INCLUDED