diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.h index 29d7d33ac4..0537068569 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.h @@ -1,101 +1,93 @@ /*=================================================================== 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 ITKDISTANCEMAPFILTER_H_ #define ITKDISTANCEMAPFILTER_H_ #include "itkImageToImageFilter.h" #include "itkImage.h" #include "mitkImage.h" namespace itk { template < class TInputImage, class TOutputImage > class DistanceMapFilter : public ImageToImageFilter { - /*! - \brief itkDistanceMapFilter - - \brief Creates a distance map from a FA skeleton image. - - \sa itkImageToImageFilter - - - \verbatim - Last contributor: $Author: vanbrugg $ - \endverbatim - */ + /** + * \brief This filter creates a map of the distance to the nearest point on a binary mask. + * + */ public: /** Typedef for input ImageType. */ typedef TInputImage InputImageType; /** Typedef for input imageType Pointer. */ typedef typename InputImageType::Pointer InputImagePointer; /** Typedef for output ImageType. */ typedef TOutputImage OutputImageType; /** Typedef for input imageType Pointer. */ typedef typename OutputImageType::Pointer OutputImagePointer; public: /** */ typedef DistanceMapFilter Self; /** Superclass */ typedef ImageToImageFilter Superclass; /** Smart Pointer */ typedef SmartPointer Pointer; /** Smart Pointer */ typedef SmartPointer ConstPointer; /** */ itkNewMacro( Self) /** Generate Data. The image will be divided into a number of pieces, a number of threads will be spawned and Threaded GenerateData() will be called in each thread. */ virtual void GenerateData(); protected: /** Constructor */ DistanceMapFilter(); /** Destructor */ virtual ~DistanceMapFilter(); protected: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkDistanceMapFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.txx index 446a3be133..f80ce05286 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkDistanceMapFilter.txx @@ -1,85 +1,79 @@ /*=================================================================== 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 _itkDistanceMapFilter_txx #define _itkDistanceMapFilter_txx #include "itkDistanceMapFilter.h" -#include "mitkProgressBar.h" #include namespace itk { template< class TInputImage, class TOutputImage > DistanceMapFilter::DistanceMapFilter() { } template< class TInputImage, class TOutputImage > DistanceMapFilter::~DistanceMapFilter() { } template< class TInputImage, class TOutputImage > void DistanceMapFilter::GenerateData() { - //----------------------------------------------------------------------// - // Progress bar // - //----------------------------------------------------------------------// - //mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 ); - typedef itk::SignedMaurerDistanceMapImageFilter DistanceFilterType; typename DistanceFilterType::Pointer dFilter = DistanceFilterType::New(); dFilter->SetInput(this->GetInput()); dFilter->SetUseImageSpacing(true); dFilter->SetSquaredDistance(false); dFilter->SetInsideIsPositive(true); dFilter->Update(); typename OutputImageType::Pointer outputImg = dFilter->GetOutput(); typedef itk::ImageRegionIterator ImageIteratorType; ImageIteratorType outIt(outputImg, outputImg->GetRequestedRegion()); outIt.GoToBegin(); while(!outIt.IsAtEnd()) { typename OutputImageType::PixelType p = outIt.Get(); p *= -1; outIt.Set(p); ++outIt; } Superclass::SetNthOutput( 0, outputImg ); } } #endif // _itkDistanceMapFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.h index f83a27f20e..163bf99444 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkProjectionFilter.h @@ -1,125 +1,155 @@ /*=================================================================== 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 ITKPROJECTIONFILTER_H_ #define ITKPROJECTIONFILTER_H_ #include "itkObject.h" #include "itkImage.h" #include "mitkImage.h" #include "mitkTbssImage.h" namespace itk { -class ProjectionFilter : public Object -{ - /*! - \brief itkProjectionFilter - - \brief Performs the Projection Step of the TBSS algorithm from Smith et al. 2006. +/** + * \brief Projection part of the TBSS pipeline + * + * This class performs the projection step of the TBSS pipeline + * (see Smith et al., 2009. http://dx.doi.org/10.1016/j.neuroimage.2006.02.024 ) + * As input it takes a binary skeleton, a distance map, and a vector image containing the image gradients + * that are typically provided by the itkSkeletonizationFilter and the itkDistanceMapFilter. */ + +class ProjectionFilter : public Object +{ + public: typedef itk::Image RealImageType; typedef itk::CovariantVector VectorType; typedef itk::Image VectorImageType; typedef itk::Image CharImageType; typedef itk::Image Float4DImageType; typedef itk::Image FloatImageType; - // typedef itk::VectorImageType SubjectVolumesType; public: /** */ typedef ProjectionFilter Self; /** Superclass */ typedef Object Superclass; /** Smart Pointer */ typedef SmartPointer Pointer; /** Smart Pointer */ typedef SmartPointer ConstPointer; /** */ itkNewMacro( Self) - /** Generate Data. The image will be divided into a number of pieces, a number of threads - will be spawned and Threaded GenerateData() will be called in each thread. */ + + /** \brief Does the actual projection */ void Project(); + + /** \brief Set the distance map + * + * Sets the distance map that decodes for every voxel the distance to the nearest point on the skeleton. + */ itkSetMacro(DistanceMap, RealImageType::Pointer) + + /** \brief Set the directions + * + * Sets the direction calculated by the TBSS skeletonization algorithm in itkSkeletonizationFilter. + */ itkSetMacro(Directions, VectorImageType::Pointer) + + /** \brief Set the binary skeleton + * + * Sets the binary skeleton that defines on which voxels must be projected. + */ itkSetMacro(Skeleton, CharImageType::Pointer) + + /** \brief Set the mask defining tubular structures on the skeleton + * + * Sets a binary mask that defines wich part of the white matter skeleton are tubular instead of sheet like. + * This is important because the a different projection method is used for sheet like structues and + * tubular structures. + */ itkSetMacro(Tube, CharImageType::Pointer) + + /** \brief Set a 4D image containing the 3D registered FA maps of all study subjects. */ itkSetMacro(AllFA, Float4DImageType::Pointer) + + /** \brief Returns a 4D image containing the skeleton projections of all subjects */ itkGetMacro(Projections, Float4DImageType::Pointer) protected: /** Constructor */ ProjectionFilter(); /** Destructor */ virtual ~ProjectionFilter(); RealImageType::Pointer m_DistanceMap; VectorImageType::Pointer m_Directions; CharImageType::Pointer m_Skeleton; CharImageType::Pointer m_Tube; Float4DImageType::Pointer m_Projections; Float4DImageType::Pointer m_AllFA; int round(float x) { if (x>0.0) return ((int) (x+0.5)); else return ((int) (x-0.5)); } protected: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkProjectionFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.h b/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.h index 26791f9095..9a46f4f167 100644 --- a/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.h +++ b/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.h @@ -1,118 +1,115 @@ /*=================================================================== 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 ITKSKELETONIZATIONFILTER_H_ #define ITKSKELETONIZATIONFILTER_H_ #include "itkImageToImageFilter.h" #include "itkImage.h" #include "mitkImage.h" #include namespace itk { -template < class TInputImage, class TOutputImage > -class SkeletonizationFilter : public ImageToImageFilter -{ - /*! - \brief itkSkeletonizationFilter - \brief Skeletonization algorithm from TBSS (Smith et al. 2006) +/** + * \brief Skeletonization part of the TBSS pipeline + * + * This class takes a 3D image (typically the mean FA image as calculated in the standard TBSS pipeline) + * and performs the non-maximum-suppression (see Smith et al., 2009. http://dx.doi.org/10.1016/j.neuroimage.2006.02.024 ) + */ - \sa itkImageToImageFilter +template < class TInputImage, class TOutputImage > +class SkeletonizationFilter : public ImageToImageFilter +{ - */ public: /** Typedef for input ImageType. */ typedef TInputImage InputImageType; typedef itk::CovariantVector VectorType; typedef itk::Image VectorImageType; /** Typedef for output ImageType. */ typedef TOutputImage OutputImageType; typedef itk::VectorImage GradientImageType; /** */ typedef SkeletonizationFilter Self; /** Superclass */ typedef ImageToImageFilter Superclass; - /** Smart Pointer */ typedef SmartPointer Pointer; - /** Smart Pointer */ typedef SmartPointer ConstPointer; - /** */ - itkNewMacro( Self); + itkNewMacro( Self) - /** Generate Data. The image will be divided into a number of pieces, a number of threads - will be spawned and Threaded GenerateData() will be called in each thread. */ + /** \brief Performs the work */ virtual void GenerateData(); + /** \brief Output the gradient image as itkVectorImage + * + * Output the gradient image by first converting it to an itk vector image + */ GradientImageType::Pointer GetGradientImage(); + /** \brief Output the gradient image as an itkImage containing vector */ VectorImageType::Pointer GetVectorImage() { return m_DirectionImage; } protected: - /** Constructor */ SkeletonizationFilter(); - /** Destructor */ virtual ~SkeletonizationFilter(); - void CalculatePerpendicularDirections(); - VectorImageType::Pointer m_DirectionImage; - int round(float x) { if (x>0.0) return ((int) (x+0.5)); else return ((int) (x-0.5)); } protected: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkSkeletonizationFilter.txx" #endif #endif diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx b/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx index 20b346e593..f42a3bcaf8 100644 --- a/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx +++ b/Modules/DiffusionImaging/Quantification/Algorithms/itkSkeletonizationFilter.txx @@ -1,313 +1,306 @@ /*=================================================================== 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 _itkSkeletonizationFilter_txx #define _itkSkeletonizationFilter_txx #include "itkSkeletonizationFilter.h" #include "mitkProgressBar.h" #include #include namespace itk { template< class TInputImage, class TOutputImage > SkeletonizationFilter::SkeletonizationFilter() { m_DirectionImage = VectorImageType::New(); } template< class TInputImage, class TOutputImage > SkeletonizationFilter::~SkeletonizationFilter() { } template< class TInputImage, class TOutputImage > void SkeletonizationFilter::GenerateData() { //----------------------------------------------------------------------// // Progress bar // //----------------------------------------------------------------------// mitk::ProgressBar::GetInstance()->AddStepsToDo( 3 ); std::cout << "Skeletonize" << std::endl; - CalculatePerpendicularDirections(); - } - - template< class TInputImage, class TOutputImage > - void SkeletonizationFilter::CalculatePerpendicularDirections() - { - const InputImageType* faImage = this->GetInput(); typename InputImageType::SizeType size = faImage->GetRequestedRegion().GetSize(); //typename RealImageType::SizeType size = m_FaImage->GetRequestedRegion().GetSize(); m_DirectionImage->SetRegions(faImage->GetRequestedRegion()); m_DirectionImage->SetDirection(faImage->GetDirection()); m_DirectionImage->SetSpacing(faImage->GetSpacing()); m_DirectionImage->SetOrigin(faImage->GetOrigin()); m_DirectionImage->Allocate(); m_DirectionImage->FillBuffer(0.0); for(int z=1; zGetPixel(ix); if(theval != 0) { /* Calculate point of gravity. We will consider each 3x3x3 neighbourhood as a unit cube. The center * point of each voxel will be a multiplicative of 1/6. The center of the unit cube is 3/6 = 1/2/ */ float cogX = 0.0; float cogY = 0.0; float cogZ = 0.0; float sum = 0.0; float l; int vecX = 0; int vecY = 0; int vecZ = 0; for(int dz=-1; dz<=1; dz++) for(int dy=-1; dy<=1; dy++) for(int dx=-1; dx<=1;dx++) { typename InputImageType::IndexType p; p[0] = x+dx; p[1] = y+dy; p[2] = z+dz; float mass = faImage->GetPixel(p); sum += mass; cogX += (float)dx*mass; cogY += (float)dy*mass; cogZ += (float)dz*mass; } cogX /= sum; cogY /= sum; cogZ /= sum; l = sqrt(cogX*cogX + cogY*cogY + cogZ*cogZ); if (l > 0.1) /* is CofG far enough away from centre voxel? */ { vecX = std::max(std::min(round(cogX/l),1),-1); vecY = std::max(std::min(round(cogY/l),1),-1); vecZ = std::max(std::min(round(cogZ/l),1),-1); } else // Find direction of max curvature { float maxcost=0, centreval=2*theval; for(int zz=0; zz<=1; zz++) // note - starts at zero as we're only searching half the voxels { for(int yy=-1; yy<=1; yy++) { for(int xx=-1; xx<=1; xx++) { if ( (zz==1) || (yy==1) || ((yy==0)&&(xx==1)) ) { float weighting = pow( (float)(xx*xx+yy*yy+zz*zz) , (float)-0.7 ); // power is arbitrary: maybe test other functions here typename InputImageType::IndexType i,j; i[0] = x+xx; i[1] = y+yy; i[2] = z+zz; j[0] = x-xx; j[1] = y-yy; j[2] = z-zz; float cost = weighting * ( centreval - (float)faImage->GetPixel(i) - (float)faImage->GetPixel(j)); if (cost>maxcost) { maxcost=cost; vecX=xx; vecY=yy; vecZ=zz; } } } } } } VectorType vec; vec[0] = vecX; vec[1] = vecY; vec[2]=vecZ; m_DirectionImage->SetPixel(ix, vec); } } mitk::ProgressBar::GetInstance()->Progress(); // Smooth m_DirectionImage and store in directionSmoothed by finding the // mode in a 3*3 neighbourhoud VectorImageType::Pointer directionSmoothed = VectorImageType::New(); directionSmoothed->SetRegions(faImage->GetRequestedRegion()); directionSmoothed->SetDirection(faImage->GetDirection()); directionSmoothed->SetSpacing(faImage->GetSpacing()); directionSmoothed->SetOrigin(faImage->GetOrigin()); directionSmoothed->Allocate(); VectorImageType::PixelType p; p[0]=0; p[1]=0; p[2]=0; directionSmoothed->FillBuffer(p); for(int z=1; zGetPixel(i); xxx = v[0]; yyy = v[1]; zzz = v[2]; localsum[(1+zzz)*9+(1+yyy)*3+1+xxx]++; localsum[(1-zzz)*9+(1-yyy)*3+1-xxx]++; } for(int zz=-1; zz<=1; zz++) for(int yy=-1; yy<=1; yy++) for(int xx=-1; xx<=1; xx++) { if (localsum[(1+zz)*9+(1+yy)*3+1+xx]>localmax) { localmax=localsum[(1+zz)*9+(1+yy)*3+1+xx]; VectorType v; v[0] = xx; v[1] = yy; v[2] = zz; directionSmoothed->SetPixel(ix, v); } } delete localsum; } m_DirectionImage = directionSmoothed; mitk::ProgressBar::GetInstance()->Progress(); // Do non-max-suppression in the direction of perp and set as output of the filter typename OutputImageType::Pointer outputImg = OutputImageType::New(); outputImg->SetRegions(faImage->GetRequestedRegion()); outputImg->SetDirection(faImage->GetDirection()); outputImg->SetSpacing(faImage->GetSpacing()); outputImg->SetOrigin(faImage->GetOrigin()); outputImg->Allocate(); outputImg->FillBuffer(0.0); for(int z=1; zGetPixel(ix); VectorType v = directionSmoothed->GetPixel(ix); typename VectorImageType::IndexType i; i[0] = x-v[0]; i[1] = y-v[1]; i[2] = z-v[2]; float min = faImage->GetPixel(i); i[0] = x+v[0]; i[1] = y+v[1]; i[2] = z+v[2]; float plus = faImage->GetPixel(i); i[0] = x-2*v[0]; i[1] = y-2*v[1]; i[2] = z-2*v[2]; float minmin = faImage->GetPixel(i); i[0] = x+2*v[0]; i[1] = y+2*v[1]; i[2] = z+2*v[2]; float plusplus = faImage->GetPixel(i); if( ((v[0]!=0) || (v[1]!=0) || (v[2]!=0)) && theval >= plus && theval > min && theval >= plusplus && theval > minmin ) { outputImg->SetPixel(ix, theval); } } Superclass::SetNthOutput( 0, outputImg ); mitk::ProgressBar::GetInstance()->Progress(); - } + // Can provide a vector image to visualize the gradient image used in the search for local maxima. template< class TInputImage, class TOutputImage > itk::VectorImage::Pointer SkeletonizationFilter::GetGradientImage() { GradientImageType::Pointer gradImg = GradientImageType::New(); if(m_DirectionImage.IsNotNull()) { gradImg->SetSpacing(m_DirectionImage->GetSpacing()); gradImg->SetOrigin(m_DirectionImage->GetOrigin()); gradImg->SetDirection(m_DirectionImage->GetDirection()); gradImg->SetRegions(m_DirectionImage->GetLargestPossibleRegion().GetSize()); gradImg->SetVectorLength(3); gradImg->Allocate(); VectorImageType::SizeType size = m_DirectionImage->GetLargestPossibleRegion().GetSize(); for(int i=0; i ix; ix[0] = i; ix[1] = j; ix[2] = k; VectorType vec = m_DirectionImage->GetPixel(ix); itk::VariableLengthVector pixel; pixel.SetSize(3); pixel.SetElement(0, vec.GetElement(0)); pixel.SetElement(1, vec.GetElement(1)); pixel.SetElement(2, vec.GetElement(2)); gradImg->SetPixel(ix, pixel); } } } } return gradImg; } } #endif // _itkSkeletonizationFilter_txx diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp index bfa25499af..7c408361cc 100644 --- a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp +++ b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.cpp @@ -1,212 +1,209 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkTractAnalyzer_cpp #define __mitkTractAnalyzer_cpp #include #include #include #include #include #include #include #include #include using namespace std; namespace mitk { - TractAnalyzer::TractAnalyzer() - { - - } + TractAnalyzer::TractAnalyzer() { } - void TractAnalyzer::BuildGraph() + void TractAnalyzer::MakeRoi() { int n = 0; if(m_PointSetNode.IsNotNull()) { n = m_PointSetNode->GetSize(); if(n==0) { QMessageBox msgBox; msgBox.setText("No points have been set yet."); msgBox.exec(); } } else{ QMessageBox msgBox; msgBox.setText("No points have been set yet."); msgBox.exec(); } std::string pathDescription = ""; std::vector< itk::Index<3> > totalPath; if(n>0) { for(int i=0; iProgress(); mitk::Point3D p = m_PointSetNode->GetPoint(i); mitk::Point3D p2 = m_PointSetNode->GetPoint(i+1); itk::Index<3> startPoint; itk::Index<3> endPoint; m_InputImage->GetGeometry()->WorldToIndex(p,startPoint); m_InputImage->GetGeometry()->WorldToIndex(p2,endPoint); MITK_INFO << "create roi"; std::vector< itk::Index<3> > path = CreateSegment(startPoint, endPoint); for(std::vector< itk::Index<3> >::iterator it = path.begin(); it != path.end(); it++) { itk::Index<3> ix = *it; if (!(ix==endPoint)) { mitk::ProgressBar::GetInstance()->Progress(); totalPath.push_back(ix); std::stringstream ss; ss << ix[0] << " " << ix[1] << " " << ix[2] << "\n"; pathDescription += ss.str(); } else { // Only when dealing with the last segment the last point should be added. This one will not occur // as the first point of the next roi segment. if(i == (n-2)) { totalPath.push_back(endPoint); std::stringstream ss; ss << endPoint[0] << " " << endPoint[1] << " " << endPoint[2] << "\n"; pathDescription += ss.str(); } } } } // save pathDescription to m_PathDescription m_PathDescription = pathDescription; FloatImageType::Pointer itkImg = FloatImageType::New(); mitk::CastToItkImage(m_InputImage, itkImg); CharImageType::Pointer roiImg = CharImageType::New(); roiImg->SetRegions(itkImg->GetLargestPossibleRegion().GetSize()); roiImg->SetOrigin(itkImg->GetOrigin()); roiImg->SetSpacing(itkImg->GetSpacing()); roiImg->SetDirection(itkImg->GetDirection()); roiImg->Allocate(); roiImg->FillBuffer(0); std::vector< itk::Index<3> > roi; std::vector< itk::Index<3> >::iterator it; for(it = totalPath.begin(); it != totalPath.end(); it++) { itk::Index<3> ix = *it; roiImg->SetPixel(ix, 1); roi.push_back(ix); } m_TbssRoi = mitk::TbssRoiImage::New(); m_TbssRoi->SetRoi(roi); m_TbssRoi->SetImage(roiImg); m_TbssRoi->InitializeFromImage(); } } std::vector< itk::Index<3> > TractAnalyzer::CreateSegment(itk::Index<3> startPoint, itk::Index<3> endPoint) { typedef itk::ShortestPathImageFilter ShortestPathFilterType; typedef itk::ShortestPathCostFunctionTbss CostFunctionType; FloatImageType::Pointer meanSkeleton; mitk::CastToItkImage(m_InputImage, meanSkeleton); // Only use the mitk image if(meanSkeleton) { CostFunctionType::Pointer costFunction = CostFunctionType::New(); costFunction->SetImage(meanSkeleton); costFunction->SetStartIndex(startPoint); costFunction->SetEndIndex(endPoint); costFunction->SetThreshold(m_Threshold); ShortestPathFilterType::Pointer pathFinder = ShortestPathFilterType::New(); pathFinder->SetCostFunction(costFunction); pathFinder->SetFullNeighborsMode(true); //pathFinder->SetCalcMode(ShortestPathFilterType::A_STAR); pathFinder->SetInput(meanSkeleton); pathFinder->SetStartIndex(startPoint); pathFinder->SetEndIndex(endPoint); pathFinder->Update(); return pathFinder->GetVectorPath(); } } } #endif diff --git a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h index 09e543dc6b..d9c3296fcf 100644 --- a/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h +++ b/Modules/DiffusionImaging/Quantification/Algorithms/mitkTractAnalyzer.h @@ -1,125 +1,152 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __mitkTractAnalyzer_h_ #define __mitkTractAnalyzer_h_ #include "QuantificationExports.h" #include -#include #include "mitkImage.h" #include "mitkImageCast.h" #include #include namespace mitk{ -/** \class TractAnalyzer - */ +/** + * \brief Creates a region of interest for tract-specific analysis of existing TBSS data + * + * This class needs a 3D image (typically a mean FA skeleton as produced by the standard TBSS pipeline of FSL) + * and a user-defined point set defining the points through which the region of interest should pass. + */ class Quantification_EXPORT TractAnalyzer { public: TractAnalyzer(); - ~TractAnalyzer() {}; + ~TractAnalyzer() {} + + /** Image type definitions */ typedef itk::Image CharImageType; typedef itk::Image FloatImageType; - typedef itk::Image ProjectionsImageType; - typedef itk::VectorImage VectorImageType; - - void BuildGraph(); - - - /* - std::vector< itk::Index<3> > GetPath() - { - return m_Path; - } - - CharImageType::Pointer GetRoiImage() - { - return m_RoiImg; - } + /** \brief Main method for region of interest calculation + * + * A region of interest is calculated adding the segments between the points on the ROI + * that was specified by the user. + */ + void MakeRoi(); -*/ + /** \brief Sets the input image + * + * The region of interest is calculated on a 3D image. This is generally the mean FA skeleton as calculated + * in the standard TBSS pipeline (see http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TBSS). + */ mitk::TbssRoiImage::Pointer GetRoiImage() { return m_TbssRoi; } + /** \brief Sets the input image + * + * The region of interest is calculated on a 3D image. This is generally the mean FA skeleton as calculated + * in the standard TBSS pipeline (see http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TBSS). + */ void SetInputImage(mitk::Image::Pointer inputImage) { m_InputImage = inputImage; } + + /** \brief Sets the user-defined point set + * + * Set the user-defined point sets. The region of interest must pass through these points. + */ void SetPointSet(mitk::PointSet::Pointer pointSet) { m_PointSetNode = pointSet; } + /** \brief Sets a lower bound for the threshold. + * + * Low fractional anisotropy values can indicate partial volume of non white matter tissue. + * This thresholds limits the search for a region of interest to voxels with a minimum value. + */ void SetThreshold(double threshold) { m_Threshold = threshold; } + + + /** \brief Returns a string with the indices of points on the region of interest + * + * The region of interest calculated by the TractAnalyzer contains a list of ITK indices. + * This method returns a string containing these indices for display in the GUI + */ std::string GetPathDescription() { return m_PathDescription; } protected: + /** \brief Calculates a segment of the region of interest + * + * The region of interest is calculated on a 3D image. This is generally the mean FA skeleton as calculated + * in the standard TBSS pipeline (see http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/TBSS). + */ std::vector< itk::Index<3> > CreateSegment(itk::Index<3> startPoint, itk::Index<3> endPoint); + /** \brief Output TbssRoiImage */ mitk::TbssRoiImage::Pointer m_TbssRoi; - //CharImageType::Pointer m_RoiImg; - + /** \brief Inputimage */ mitk::Image::Pointer m_InputImage; + /** \brief Threshold for ROI search */ double m_Threshold; - //std::vector< itk::Index<3> > m_Path; - + /** \brief User defined point set */ mitk::PointSet::Pointer m_PointSetNode; + /** \brief Path description in as string for display in GUI */ std::string m_PathDescription; private: }; } #endif //__itkTractAnalyzer_h_ diff --git a/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.cpp b/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.cpp index f184dac8be..c4c3b939b5 100644 --- a/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.cpp +++ b/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.cpp @@ -1,132 +1,126 @@ /*=================================================================== 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 __mitkTbssImporter_cpp #define __mitkTbssImporter_cpp #include "mitkTbssImporter.h" -#include -#include -#include "itkNrrdImageIO.h" -#include "mitkImageAccessByItk.h" namespace mitk { mitk::TbssImage::Pointer mitk::TbssImporter::Import() { // read all images with all_*.nii.gz mitk::TbssImage::Pointer tbssImg = mitk::TbssImage::New(); m_Data = DataImageType::New(); mitk::Geometry3D* geo = m_InputVolume->GetGeometry(); mitk::Vector3D spacing = geo->GetSpacing(); mitk::Point3D origin = geo->GetOrigin(); //Size size DataImageType::SizeType dataSize; dataSize[0] = m_InputVolume->GetDimension(0); dataSize[1] = m_InputVolume->GetDimension(1); dataSize[2] = m_InputVolume->GetDimension(2); m_Data->SetRegions(dataSize); // Set spacing DataImageType::SpacingType dataSpacing; dataSpacing[0] = spacing[0]; dataSpacing[1] = spacing[1]; dataSpacing[2] = spacing[2]; m_Data->SetSpacing(dataSpacing); DataImageType::PointType dataOrigin; dataOrigin[0] = origin[0]; dataOrigin[1] = origin[1]; dataOrigin[2] = origin[2]; m_Data->SetOrigin(dataOrigin); //Direction must be set DataImageType::DirectionType dir; const itk::Transform* transform3D = geo->GetParametricTransform(); itk::Transform::ParametersType p = transform3D->GetParameters(); int t=0; for(int i=0; i<3; i++) { for(int j=0; j<3; j++) { dir[i][j] = p[t]; // row-major order (where the column index varies the fastest) t++; } } m_Data->SetDirection(dir); // Set the length to one because otherwise allocate fails. Should be changed when groups/measurements are added m_Data->SetVectorLength(m_InputVolume->GetDimension(3)); m_Data->Allocate(); for(int i=0; i ix; ix[0] = i; ix[1] = j; ix[2] = k; itk::VariableLengthVector pixel = m_Data->GetPixel(ix); for(int z=0; zGetPixelValueByIndex(ix, z); pixel.SetElement(z, value); } m_Data->SetPixel(ix, pixel); } } } -// mitk::CastToTbssImage(m_Data.GetPointer(), tbssImg); - tbssImg->SetGroupInfo(m_Groups); tbssImg->SetMeasurementInfo(m_MeasurementInfo); tbssImg->SetImage(m_Data); tbssImg->InitializeFromVectorImage(); return tbssImg; } } #endif // __mitkTbssImporter_cpp diff --git a/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.h b/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.h index 4fb4254602..c5b6390ad2 100644 --- a/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.h +++ b/Modules/DiffusionImaging/Quantification/IODataStructures/TbssImages/mitkTbssImporter.h @@ -1,104 +1,96 @@ /*=================================================================== 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 __mitkTbssImporter_h #define __mitkTbssImporter_h #include "mitkCommon.h" #include "mitkFileReader.h" #include "itkImage.h" #include "itkVectorImage.h" -#include "itkImageFileReader.h" + #include "mitkTbssImage.h" #include "QuantificationExports.h" - namespace mitk { - //template - class Quantification_EXPORT TbssImporter : public itk::Object { - public: - // typedef TPixelType PixelType; - typedef itk::VectorImage DataImageType; // type of the 3d vector image containing the skeletonized images - typedef itk::VectorImage VectorImageType; // Datatype of the tbss gradient images - typedef itk::Image FloatImage4DType; - typedef itk::ImageFileReader FileReaderType4D; - typedef itk::ImageFileReader VectorReaderType; +/** + * \brief Converts FSL TBSS data (4D skeleton projection images) to a NRRD image with meta data. + * + * The TBSS pipeline of FSL produces a 4D image containing the 3D skeleton projections of all individuals. + * This class converts the FSL Nifty image to NRRD and adds information about the type of measurement and the study groups. + */ - typedef itk::Image FloatImage3DType; - typedef itk::ImageFileReader FileReaderType3D; + class Quantification_EXPORT TbssImporter : public itk::Object { + public: + // type of the 3d vector image containing the skeletonized images + typedef itk::VectorImage DataImageType; mitkClassMacro( TbssImporter, Object ) itkNewMacro(Self) - - + /* \brief Converts the FSL Nifty to NRRD and adds the meta data */ mitk::TbssImage::Pointer Import(); + + /* \brief Group info is set by providing a vector with pairs of group name and number*/ void SetGroupInfo(std::vector< std::pair > groups) { m_Groups = groups; } - std::vector< std::pair > GetGroupInfo() - { - return m_Groups; - } + /* \brief Used to indicate the type of measurement */ void SetMeasurementInfo(std::string s) { m_MeasurementInfo = s; } - std::string GetMeasurementInfo() - { - return m_MeasurementInfo; - } - + /* \brief Sets the FSL import volume */ void SetImportVolume(mitk::Image::Pointer inputVolume) { m_InputVolume = inputVolume; } protected: TbssImporter(){} virtual ~TbssImporter(){} DataImageType::Pointer m_Data; + std::vector< std::pair > m_Groups; std::string m_MeasurementInfo; - mitk::Image::Pointer m_InputVolume; }; } #endif // __mitkTbssImporter_h diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp index f9007d45d2..118e7e599c 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.cpp @@ -1,973 +1,968 @@ /*=================================================================== 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(); mitk::PlaneGeometry* startGeometry2D = dynamic_cast( const_cast(startRoi->GetGeometry2D()) ); 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(); 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(NULL); lines->GetNextCell ( numPointsInCell, pointsInCell ); int startId = 0; int endId = 0; float minDistStart = std::numeric_limits::max(); float minDistEnd = std::numeric_limits::max(); 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*/ 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 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]; 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]; } 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]; 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 *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 double *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 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]; 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]; } 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]; 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 *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 double *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); } //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 PlotFiberBundles(resampledTracts, img, avg); m_CurrentTracts = resampledTracts; } 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) return; DoPlotFiberBundles(m_Fib, m_CurrentImage, m_CurrentStartRoi, m_CurrentEndRoi, avg, 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) { MITK_INFO << "make profiles!"; std::vector< std::vector > profiles; //No results were preprocessed, so they must be calculated now. if(preprocessed == "") { // Iterate through the 4th dim (corresponding to subjects) // and create a profile for every subject int size = m_Projections->GetVectorLength(); for(int s=0; s profile; RoiType::iterator it; it = m_Roi.begin(); while(it != m_Roi.end()) { itk::Index<3> ix = *it; profile.push_back(m_Projections->GetPixel(ix).GetElement(s)); it++; } int pSize = profile.size(); profiles.push_back(profile); } } else{ // Use preprocessed results std::ifstream file(preprocessed.c_str()); if(file.is_open()) { std::string line; while(getline(file,line)) { std::vector tokens; Tokenize(line, tokens); std::vector::iterator it; it = tokens.begin(); std::vector< double > profile; while(it != tokens.end()) { std::string s = *it; profile.push_back (atof( s.c_str() ) ); ++it; } profiles.push_back(profile); } } } 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(); m_Vals.clear(); std::vector v1; std::vector > groupProfiles = CalculateGroupProfiles(preprocessed); std::vector xAxis; for(int i=0; iSetPlotTitle( title.c_str() ); QPen pen( Qt::SolidLine ); pen.setWidth(2); std::vector< std::pair >::iterator it; it = m_Groups.begin(); int c = 0; //the current profile number QColor colors[4] = {Qt::green, Qt::blue, Qt::yellow, Qt::red}; while(it != m_Groups.end() && groupProfiles.size() > 0) { std::pair< std::string, int > group = *it; pen.setColor(colors[c]); int curveId = this->InsertCurve( group.first.c_str() ); this->SetCurveData( curveId, xAxis, groupProfiles.at(c) ); this->SetCurvePen( curveId, pen ); c++; it++; } QwtLegend *legend = new QwtLegend; this->SetLegend(legend, QwtPlot::RightLegend, 0.5); std::cout << m_Measure << std::endl; this->m_Plot->setAxisTitle(0, m_Measure.c_str()); this->m_Plot->setAxisTitle(3, "Position"); this->Replot(); } void QmitkTbssRoiAnalysisWidget::PlotFiberBundles(TractContainerType tracts, mitk::Image *img, bool avg) { this->Clear(); std::vector::iterator it = tracts.begin(); // Match points on tracts. Take the smallest tract and match all others on this one /* int min = std::numeric_limits::max(); TractType smallestTract; while(it != tracts.end()) { TractType tract = *it; if(tract.size() correspondingIndices; TractType correspondingPoints; for(int i=0; i::max(); int correspondingIndex = 0; PointType correspondingPoint; // Search for the point on the second tract with the smallest distance // to p and memorize it for(int j=0; j > 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() ); 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() ); 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::Boxplots() -{ - this->Clear(); -} - void QmitkTbssRoiAnalysisWidget::drawBar(int x) { m_Plot->detachItems(QwtPlotItem::Rtti_PlotMarker, true); QwtPlotMarker *mX = new QwtPlotMarker(); //mX->setLabel(QString::fromLatin1("selected point")); mX->setLabelAlignment(Qt::AlignLeft | Qt::AlignBottom); mX->setLabelOrientation(Qt::Vertical); mX->setLineStyle(QwtPlotMarker::VLine); mX->setLinePen(QPen(Qt::black, 0, Qt::SolidLine)); mX->setXValue(x); mX->attach(m_Plot); this->Replot(); } QmitkTbssRoiAnalysisWidget::~QmitkTbssRoiAnalysisWidget() { delete m_PlotPicker; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h index 3002f61d35..6c85b4f9d7 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/QmitkTbssRoiAnalysisWidget.h @@ -1,223 +1,217 @@ /*=================================================================== 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 "QmitkHistogram.h" #include "QmitkExtExports.h" -#include "mitkImage.h" #include "mitkPlanarFigure.h" #include "itkVectorImage.h" #include #include -//#include - -//#include -#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 Widget for displaying boxplots - * framework + * \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 PlotFiberBundles(TractContainerType tracts, mitk::Image* img, bool avg=false); - void Boxplots(); - + /* \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; } - QwtPlot* GetPlot() - { - return m_Plot; - } - - QwtPlotPicker* m_PlotPicker; - + /* \brief Draws a bar to indicate were the user clicked in the plot */ + void drawBar(int x); - 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); 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::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); }; #endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.h index b13712592e..739891fd75 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssMetaTableModel.h @@ -1,54 +1,54 @@ /*=================================================================== 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 QmitkTbssMetaTableModel_h #define QmitkTbssMetaTableModel_h #include #include #include #include #include class QmitkTbssMetaTableModel : public QAbstractTableModel { - //Q_OBJECT + // TableModel to model pairs of group description and group number public: QmitkTbssMetaTableModel(QObject *parent=0); QmitkTbssMetaTableModel(QList< QPair > listofPairs, QObject *parent=0); int rowCount(const QModelIndex &parent) const; int columnCount(const QModelIndex &parent) const; QVariant data(const QModelIndex &index, int role) const; QVariant headerData(int section, Qt::Orientation orientation, int role) const; Qt::ItemFlags flags(const QModelIndex &index) const; bool setData(const QModelIndex &index, const QVariant &value, int role=Qt::EditRole); bool insertRows(int position, int rows, const QModelIndex &index=QModelIndex()); bool removeRows(int position, int rows, const QModelIndex &index=QModelIndex()); QList< QPair > getList(); private: QList< QPair > listOfPairs; }; #endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.cpp index 5a69a8175f..755023f702 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.cpp @@ -1,453 +1,435 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ -// Blueberry -#include -#include -#include -#include -#include - // Qmitk #include "QmitkTbssSkeletonizationView.h" -#include -#include -#include #include #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() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) { } QmitkTbssSkeletonizationView::~QmitkTbssSkeletonizationView() { } void QmitkTbssSkeletonizationView::OnSelectionChanged(std::vector nodes) { //datamanager selection changed if (!this->IsActivated()) return; bool found3dImage = false; bool found4dImage = false; // iterate selection for ( int i=0; iGetData(); if(nodeData) { if(QString("Image").compare(nodeData->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodeData); if(img->GetDimension() == 3) { found3dImage = true; } else if(img->GetDimension() == 4) { found4dImage = 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(); } - m_IsInitialized = false; - - } void QmitkTbssSkeletonizationView::Activated() { QmitkFunctionality::Activated(); - - berry::ISelection::ConstPointer sel( - this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); - //m_CurrentSelection = sel.Cast(); - //m_SelListener.Cast()->DoSelectionChanged(sel); } void QmitkTbssSkeletonizationView::Deactivated() { QmitkFunctionality::Deactivated(); } 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::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkTbssSkeletonizationView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkTbssSkeletonizationView::Skeletonize() { typedef itk::SkeletonizationFilter SkeletonisationFilterType; SkeletonisationFilterType::Pointer skeletonizer = SkeletonisationFilterType::New(); std::vector nodes = this->GetDataManagerSelection(); mitk::Image::Pointer meanImage = mitk::Image::New(); std::string name = ""; for ( int i=0; iGetData(); if(nodeData) { if(QString("Image").compare(nodeData->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodeData); if(img->GetDimension() == 3) { meanImage = img; name = nodes[i]->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(); std::vector nodes = this->GetDataManagerSelection(); mitk::Image::Pointer meanImage = mitk::Image::New(); mitk::Image::Pointer subjects = mitk::Image::New(); for ( int i=0; iGetData(); if(nodeData) { if(QString("Image").compare(nodeData->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodeData); if(img->GetDimension() == 3) { meanImage = img; } else if(img->GetDimension() == 4) { subjects = img; } } } } Float4DImageType::Pointer allFA = ConvertToItk(subjects); // 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); } typedef itk::ImageFileReader< CharImageType > CharReaderType; CharReaderType::Pointer reader = CharReaderType::New(); reader->SetFileName("/local/testing/LowerCingulum_1mm.nii.gz"); reader->Update(); CharImageType::Pointer cingulum = reader->GetOutput(); ProjectionFilterType::Pointer projectionFilter = ProjectionFilterType::New(); projectionFilter->SetDistanceMap(distanceMap); projectionFilter->SetDirections(directionImg); projectionFilter->SetAllFA(allFA); projectionFilter->SetTube(cingulum); 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 GetDefaultDataStorage()->Add( result ); } Float4DImageType::Pointer QmitkTbssSkeletonizationView::ConvertToItk(mitk::Image::Pointer image) { Float4DImageType::Pointer output = Float4DImageType::New(); mitk::Geometry3D* 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); // iterate through the subjects and copy data to output for(int t=0; tGetDimension(0); x++) { for(int y=0; yGetDimension(1); y++) { for(int z=0; zGetDimension(2); z++) { itk::Index<4> ix4; ix4[0] = x; ix4[1] = y; ix4[2] = z; ix4[3] = t; mitk::Index3D ix; ix[0] = x; ix[1] = y; ix[2] = z; output->SetPixel(ix4, image->GetPixelValueByIndex(ix, t)); } } } } } return output; } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.h index 2c461bc102..683c7bdd1b 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTbssSkeletonizationView.h @@ -1,115 +1,95 @@ /*=================================================================== 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 QmitkTbssSkeletonizationView_h #define QmitkTbssSkeletonizationView_h -// itk -#include "itkImage.h" -#include "itkVectorImage.h" - #include - #include "ui_QmitkTbssSkeletonizationViewControls.h" typedef itk::Image FloatImageType; typedef itk::Image CharImageType; typedef itk::Image Float4DImageType; typedef itk::CovariantVector VectorType; typedef itk::Image DirectionImageType; -struct TbssSkeletonizationSelListener; - - /*! - \brief QmitkTbssSkeletonizationView - - \warning This application module is not yet documented. Use "svn blame/praise/annotate" and ask the author to provide basic documentation. - - \sa QmitkFunctionalitymitkTbssWorkspaceManager - \ingroup Functionalities + * \brief Implementation of the core functionality of TBSS. + * This plugin provides the core functionality of TBSS (see Smith et al., 2009. http://dx.doi.org/10.1016/j.neuroimage.2006.02.024) + * It can skeletonize a mean FA image and calculate the projection of all individual subjects to this skeleton. */ + class QmitkTbssSkeletonizationView : public QmitkFunctionality { - - friend struct TbssSkeletonizationSelListener; - - - - // this is needed for all Qt objesetupUicts that should have a Qt meta-object - // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkTbssSkeletonizationView(); virtual ~QmitkTbssSkeletonizationView(); virtual void CreateQtPartControl(QWidget *parent); - /// \brief Creation of the connections of main and control widget + //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: + /* \brief Perform skeletonization only */ void Skeletonize(); + + // Perform skeletonization and Projection of subject data to the skeleton void Project(); protected: - /// \brief called by QmitkFunctionality when DataManager's selection has changed + //brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); - void InitPointsets(); - - void SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name); - - bool m_IsInitialized; - Ui::QmitkTbssSkeletonizationViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; void AddToDataStorage(mitk::Image* img, std::string name); Float4DImageType::Pointer ConvertToItk(mitk::Image::Pointer image); }; #endif // _QMITKTbssSkeletonizationVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp index 79bc0e16dc..db6de5555f 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.cpp @@ -1,1315 +1,1132 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ -// Blueberry -#include "berryIWorkbenchWindow.h" -#include "berryIWorkbenchPage.h" -#include "berryISelectionService.h" -#include "berryConstants.h" -#include "berryPlatformUI.h" - // Qmitk #include "QmitkTractbasedSpatialStatisticsView.h" #include "QmitkStdMultiWidget.h" #include "mitkDataNodeObject.h" #include // Qt #include -#include -#include #include -#include -#include -#include -#include #include -#include -#include - -#include -#include -#include - -#include #include #include -#include #include #include -#include -#include "vtkFloatArray.h" -#include "vtkLinearTransform.h" #include "vtkPoints.h" -#include "mitkSurface.h" -#include -#include "vtkArrowSource.h" -#include "vtkUnstructuredGrid.h" -#include "vtkPointData.h" #include #include -#include -#include -#include - -#include - -#include "mitkITKImageImport.h" -// #include "mitkImageMapperGL2D.h" -#include "mitkVolumeDataVtkMapper3D.h" -#include "mitkImageAccessByItk.h" -#include "mitkTensorImage.h" - -#include "itkDiffusionTensor3D.h" - - -#define SEARCHSIGMA 10 // length in linear voxel dimens - -#define MAXSEARCHLENGTH (3*SEARCHSIGMA) const std::string QmitkTractbasedSpatialStatisticsView::VIEW_ID = "org.mitk.views.tractbasedspatialstatistics"; using namespace berry; QmitkTractbasedSpatialStatisticsView::QmitkTractbasedSpatialStatisticsView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) { } QmitkTractbasedSpatialStatisticsView::~QmitkTractbasedSpatialStatisticsView() { } void QmitkTractbasedSpatialStatisticsView::PerformChange() { m_Controls->m_RoiPlotWidget->ModifyPlot(m_Controls->m_Segments->value(), m_Controls->m_Average->isChecked()); } void QmitkTractbasedSpatialStatisticsView::OnSelectionChanged(std::vector nodes) { //datamanager selection changed if (!this->IsActivated()) return; - +/* // Get DataManagerSelection if (!this->GetDataManagerSelection().empty()) { mitk::DataNode::Pointer sourceImageNode = this->GetDataManagerSelection().front(); mitk::Image::Pointer sourceImage = dynamic_cast(sourceImageNode->GetData()); if (!sourceImage) { m_Controls->m_TbssImageLabel->setText( QString( sourceImageNode->GetName().c_str() ) + " is no image" ); return; } // set Text m_Controls->m_TbssImageLabel->setText( QString( sourceImageNode->GetName().c_str() ) + " (" + QString::number(sourceImage->GetDimension()) + "D)" ); } else { m_Controls->m_TbssImageLabel->setText("Please select an image"); } +*/ + // 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) { this->PlotFiberBundle(fib, img, start, end); } else if(found3dImage == true && foundFiberBundle) { this->PlotFiberBundle(fib, img); } if(found3dImage) { this->InitPointsets(); } this->m_Controls->m_Cut->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi); this->m_Controls->m_SegmentLabel->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); this->m_Controls->m_Segments->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); this->m_Controls->m_Average->setEnabled(foundFiberBundle && foundStartRoi && foundEndRoi && found3dImage); } void QmitkTractbasedSpatialStatisticsView::InitPointsets() { // Check if PointSetStart exsits, if not create it. m_P1 = this->GetDefaultDataStorage()->GetNamedNode("PointSetNode"); if (m_PointSetNode) { //m_PointSetNode = dynamic_cast(m_P1->GetData()); return; } if ((!m_P1) || (!m_PointSetNode)) { // create new ones m_PointSetNode = mitk::PointSet::New(); m_P1 = mitk::DataNode::New(); m_P1->SetData( m_PointSetNode ); m_P1->SetProperty( "name", mitk::StringProperty::New( "PointSet" ) ); m_P1->SetProperty( "opacity", mitk::FloatProperty::New( 1 ) ); m_P1->SetProperty( "helper object", mitk::BoolProperty::New(true) ); // CHANGE if wanted m_P1->SetProperty( "pointsize", mitk::FloatProperty::New( 0.1 ) ); m_P1->SetColor( 1.0, 0.0, 0.0 ); this->GetDefaultDataStorage()->Add(m_P1); m_Controls->m_PointWidget->SetPointSetNode(m_P1); m_Controls->m_PointWidget->SetMultiWidget(GetActiveStdMultiWidget()); } } void QmitkTractbasedSpatialStatisticsView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkTractbasedSpatialStatisticsViewControls; m_Controls->setupUi( parent ); this->CreateConnections(); } - m_IsInitialized = false; - // Table for the FSL TBSS import m_GroupModel = new QmitkTbssTableModel(); m_Controls->m_GroupInfo->setModel(m_GroupModel); } void QmitkTractbasedSpatialStatisticsView::Activated() { QmitkFunctionality::Activated(); - berry::ISelection::ConstPointer sel( - this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); - //m_CurrentSelection = sel.Cast(); - //m_SelListener.Cast()->DoSelectionChanged(sel); } void QmitkTractbasedSpatialStatisticsView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkTractbasedSpatialStatisticsView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_CreateRoi), SIGNAL(clicked()), this, SLOT(CreateRoi()) ); connect( (QObject*)(m_Controls->m_ImportFsl), SIGNAL(clicked()), this, SLOT(TbssImport()) ); connect( (QObject*)(m_Controls->m_AddGroup), SIGNAL(clicked()), this, SLOT(AddGroup()) ); connect( (QObject*)(m_Controls->m_RemoveGroup), SIGNAL(clicked()), this, SLOT(RemoveGroup()) ); connect( (QObject*)(m_Controls->m_Clipboard), SIGNAL(clicked()), this, SLOT(CopyToClipboard()) ); connect( m_Controls->m_RoiPlotWidget->m_PlotPicker, SIGNAL(selected(const QwtDoublePoint&)), SLOT(Clicked(const QwtDoublePoint&) ) ); connect( m_Controls->m_RoiPlotWidget->m_PlotPicker, SIGNAL(moved(const QwtDoublePoint&)), SLOT(Clicked(const QwtDoublePoint&) ) ); connect( (QObject*)(m_Controls->m_Cut), SIGNAL(clicked()), this, SLOT(Cut()) ); connect( (QObject*)(m_Controls->m_Average), SIGNAL(stateChanged(int)), this, SLOT(PerformChange()) ); connect( (QObject*)(m_Controls->m_Segments), SIGNAL(valueChanged(int)), this, SLOT(PerformChange()) ); } } void QmitkTractbasedSpatialStatisticsView::CopyToClipboard() { if(m_Controls->m_RoiPlotWidget->IsPlottingFiber()) { // Working with fiber bundles std::vector > profiles = m_Controls->m_RoiPlotWidget->GetIndividualProfiles(); - - QString clipboardText; for (std::vector >::iterator it = profiles.begin(); it != profiles.end(); ++it) { for (std::vector::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); } clipboardText.append(QString("\n")); } if(m_Controls->m_Average->isChecked()) { std::vector averages = m_Controls->m_RoiPlotWidget->GetAverageProfile(); clipboardText.append(QString("\nAverage\n")); for (std::vector::iterator it2 = averages.begin(); it2 != averages.end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); } } QApplication::clipboard()->setText(clipboardText, QClipboard::Clipboard); } else{ // Working with TBSS Data if(m_Controls->m_Average->isChecked()) { std::vector > vals = m_Controls->m_RoiPlotWidget->GetVals(); QString clipboardText; for (std::vector >::iterator it = vals.begin(); it != vals.end(); ++it) { for (std::vector::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); double d = *it2; std::cout << d <setText(clipboardText, QClipboard::Clipboard); } else { std::vector > vals = m_Controls->m_RoiPlotWidget->GetIndividualProfiles(); QString clipboardText; for (std::vector >::iterator it = vals.begin(); it != vals.end(); ++it) { for (std::vector::iterator it2 = (*it).begin(); it2 != (*it).end(); ++it2) { clipboardText.append(QString("%1 \t").arg(*it2)); double d = *it2; std::cout << d <setText(clipboardText, QClipboard::Clipboard); } } } void QmitkTractbasedSpatialStatisticsView::RemoveGroup() { QTableView *temp = static_cast(m_Controls->m_GroupInfo); - // QSortFilterProxyModel *proxy = static_cast(temp->model()); + QItemSelectionModel *selectionModel = temp->selectionModel(); QModelIndexList indices = selectionModel->selectedRows(); QModelIndex index; foreach(index, indices) { int row = index.row(); m_GroupModel->removeRows(row, 1, QModelIndex()); } } -std::string QmitkTractbasedSpatialStatisticsView::ReadFile(std::string whatfile) -{ - std::string s = "Select a" + whatfile; - QFileDialog* w = new QFileDialog(this->m_Controls->m_ImportFsl, QString(s.c_str()) ); - w->setFileMode(QFileDialog::ExistingFiles); - w->setDirectory("/home"); - - if(whatfile == "gradient image") - { - w->setNameFilter("Tbss gradient images (*.tgi)"); - } - // RETRIEVE SELECTION - if ( w->exec() != QDialog::Accepted ) - { - return ""; - MITK_INFO << "Failed to load"; - } - - QStringList filenames = w->selectedFiles(); - if (filenames.size() > 0) - { - std::string retval = filenames.at(0).toStdString(); - return retval; - } - return ""; -} void QmitkTractbasedSpatialStatisticsView::AddGroup() { QString group("Group"); int number = 0; QPair pair(group, number); QList< QPair >list = m_GroupModel->getList(); if(!list.contains(pair)) { m_GroupModel->insertRows(0, 1, QModelIndex()); QModelIndex index = m_GroupModel->index(0, 0, QModelIndex()); m_GroupModel->setData(index, group, Qt::EditRole); index = m_GroupModel->index(0, 1, QModelIndex()); m_GroupModel->setData(index, number, Qt::EditRole); } - else - { - //QMessageBox::information(this, "Duplicate name"); - } - - } void QmitkTractbasedSpatialStatisticsView::TbssImport() { // Read groups from the interface mitk::TbssImporter::Pointer importer = mitk::TbssImporter::New(); QList< QPair >list = m_GroupModel->getList(); + if(list.size() == 0) { QMessageBox msgBox; msgBox.setText("No study group information has been set yet."); msgBox.exec(); return; } std::vector < std::pair > groups; for(int i=0; i pair = list.at(i); std::string s = pair.first.toStdString(); int n = pair.second; std::pair p; p.first = s; p.second = n; groups.push_back(p); } importer->SetGroupInfo(groups); std::string minfo = m_Controls->m_MeasurementInfo->text().toStdString(); importer->SetMeasurementInfo(minfo); std::string name = ""; std::vector nodes = this->GetDataManagerSelection(); for ( int i=0; iGetData()->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodes[i]->GetData()); if(img->GetDimension() == 4) { importer->SetImportVolume(img); name = nodes[i]->GetName(); } } } mitk::TbssImage::Pointer tbssImage; tbssImage = importer->Import(); name += "_tbss"; AddTbssToDataStorage(tbssImage, name); } void QmitkTractbasedSpatialStatisticsView::AddTbssToDataStorage(mitk::Image* image, std::string name) { mitk::LevelWindow levelwindow; levelwindow.SetAuto( image ); mitk::LevelWindowProperty::Pointer levWinProp = mitk::LevelWindowProperty::New(); levWinProp->SetLevelWindow( levelwindow ); mitk::DataNode::Pointer result = mitk::DataNode::New(); result->SetProperty( "name", mitk::StringProperty::New(name) ); result->SetData( image ); result->SetProperty( "levelwindow", levWinProp ); // add new image to data storage and set as active to ease further processing GetDefaultDataStorage()->Add( result ); // show the results mitk::RenderingManager::GetInstance()->RequestUpdateAll(); } + void QmitkTractbasedSpatialStatisticsView::Clicked(const QwtDoublePoint& pos) { int index = (int)pos.x(); if(m_Roi.size() > 0 && m_CurrentGeometry != NULL && !m_Controls->m_RoiPlotWidget->IsPlottingFiber() ) { index = std::min( (int)m_Roi.size()-1, std::max(0, index) ); itk::Index<3> ix = m_Roi.at(index); mitk::Vector3D i; i[0] = ix[0]; i[1] = ix[1]; i[2] = ix[2]; mitk::Vector3D w; m_CurrentGeometry->IndexToWorld(i, w); mitk::Point3D origin = m_CurrentGeometry->GetOrigin(); mitk::Point3D p; p[0] = w[0] + origin[0]; p[1] = w[1] + origin[1]; p[2] = w[2] + origin[2]; m_MultiWidget->MoveCrossToPosition(p); m_Controls->m_RoiPlotWidget->drawBar(index); } else if(m_Controls->m_RoiPlotWidget->IsPlottingFiber() ) { mitk::Point3D point = m_Controls->m_RoiPlotWidget->GetPositionInWorld(index); m_MultiWidget->MoveCrossToPosition(point); } } void QmitkTractbasedSpatialStatisticsView::Cut() { mitk::BaseData* fibData = m_CurrentFiberNode->GetData(); mitk::FiberBundleX* fib = static_cast(fibData); mitk::BaseData* startData = m_CurrentStartRoi->GetData(); mitk::PlanarFigure* startRoi = static_cast(startData); mitk::PlaneGeometry* startGeometry2D = dynamic_cast( const_cast(startRoi->GetGeometry2D()) ); mitk::BaseData* endData = m_CurrentEndRoi->GetData(); mitk::PlanarFigure* endRoi = static_cast(endData); mitk::PlaneGeometry* endGeometry2D = dynamic_cast( const_cast(endRoi->GetGeometry2D()) ); mitk::Point3D startCenter = startRoi->GetWorldControlPoint(0); //center Point of start roi mitk::Point3D endCenter = endRoi->GetWorldControlPoint(0); //center Point of end roi mitk::FiberBundleX::Pointer inStart = fib->ExtractFiberSubset(startRoi); mitk::FiberBundleX::Pointer inBoth = inStart->ExtractFiberSubset(endRoi); int num = inBoth->GetNumFibers(); vtkSmartPointer fiberPolyData = inBoth->GetFiberPolyData(); vtkCellArray* lines = fiberPolyData->GetLines(); lines->InitTraversal(); // initialize new vtk polydata vtkSmartPointer points = vtkSmartPointer::New(); vtkSmartPointer polyData = vtkSmartPointer::New(); vtkSmartPointer cells = vtkSmartPointer::New(); int pointIndex=0; // find start and endpoint for( int fiberID( 0 ); fiberID < num; fiberID++ ) { vtkIdType numPointsInCell(0); vtkIdType* pointsInCell(NULL); lines->GetNextCell ( numPointsInCell, pointsInCell ); int startId = 0; int endId = 0; float minDistStart = std::numeric_limits::max(); float minDistEnd = std::numeric_limits::max(); vtkSmartPointer polyLine = vtkSmartPointer::New(); int lineIndex=0; for( int pointInCellID( 0 ); pointInCellID < numPointsInCell ; pointInCellID++) { double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); mitk::Point3D point; point[0] = p[0]; point[1] = p[1]; point[2] = p[2]; float distanceToStart = point.EuclideanDistanceTo(startCenter); float distanceToEnd = point.EuclideanDistanceTo(endCenter); if(distanceToStart < minDistStart) { minDistStart = distanceToStart; startId = pointInCellID; } if(distanceToEnd < minDistEnd) { minDistEnd = distanceToEnd; endId = pointInCellID; } } /* We found the start and end points of of the part that should be plottet for the current fiber. now we need to plot them. If the endId is smaller than the startId the plot order must be reversed*/ if(startId < endId) { double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId+1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool secondOnPositive = startGeometry2D->IsAbove(pSecond); // Calculate intersection with the plane mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId+1 lie on the same side of the plane, so we need need startId-1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; double insertPoint[3]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; // First insert the intersection with the start roi points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; if(! (secondOnPositive ^ startOnPositive) ) { /* StartId and startId+1 lie on the same side of the plane so startId is also part of the ROI*/ double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); points->InsertNextPoint(start); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Insert the rest up and to including endId-1 for( int pointInCellID( startId+1 ); pointInCellID < endId ; pointInCellID++) { // create new polyline for new polydata double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); points->InsertNextPoint(p); // add point to line polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } /* endId must be included if endId and endId-1 lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; bool lastOnPositive = endGeometry2D->IsAbove(pLast); bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); normal = endGeometry2D->GetNormal(); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId-1 lie on the same side of the plane, so we need need endId+1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure endId is part of the roi so can also be included here*/ p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; double *end = fiberPolyData->GetPoint( pointsInCell[endId] ); points->InsertNextPoint(end); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; //Insert the Last Point (intersection with the end roi) points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Need to reverse walking order else{ double *p = fiberPolyData->GetPoint( pointsInCell[ startId ] ); mitk::Vector3D p0; p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ startId-1 ] ); mitk::Vector3D p1; p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; // Check if p and p2 are both on the same side of the plane mitk::Vector3D normal = startGeometry2D->GetNormal(); mitk::Point3D pStart; pStart[0] = p0[0]; pStart[1] = p0[1]; pStart[2] = p0[2]; bool startOnPositive = startGeometry2D->IsAbove(pStart); mitk::Point3D pSecond; pSecond[0] = p1[0]; pSecond[1] = p1[1]; pSecond[2] = p1[2]; bool secondOnPositive = startGeometry2D->IsAbove(pSecond); // Calculate intersection with the plane mitk::Vector3D onPlane; onPlane[0] = startCenter[0]; onPlane[1] = startCenter[1]; onPlane[2] = startCenter[2]; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId-1 lie on the same side of the plane, so we need need startId+1 to calculate the intersection with the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ startId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; } double d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); mitk::Vector3D newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; double insertPoint[3]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; // First insert the intersection with the start roi points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; if(! (secondOnPositive ^ startOnPositive) ) { /* startId and startId-1 lie on the same side of the plane so endId is also part of the ROI*/ double *start = fiberPolyData->GetPoint( pointsInCell[startId] ); points->InsertNextPoint(start); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // Insert the rest up and to including endId-1 for( int pointInCellID( startId-1 ); pointInCellID > endId ; pointInCellID--) { // create new polyline for new polydata double *p = fiberPolyData->GetPoint( pointsInCell[ pointInCellID ] ); points->InsertNextPoint(p); // add point to line polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } /* startId must be included if startId and startId+ lie on the same side of the plane defined by endRoi*/ p = fiberPolyData->GetPoint( pointsInCell[ endId ] ); p0[0] = p[0]; p0[1] = p[1]; p0[2] = p[2]; p = fiberPolyData->GetPoint( pointsInCell[ endId+1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; mitk::Point3D pLast; pLast[0] = p0[0]; pLast[1] = p0[1]; pLast[2] = p0[2]; bool lastOnPositive = endGeometry2D->IsAbove(pLast); mitk::Point3D pBeforeLast; pBeforeLast[0] = p1[0]; pBeforeLast[1] = p1[1]; pBeforeLast[2] = p1[2]; bool secondLastOnPositive = endGeometry2D->IsAbove(pBeforeLast); onPlane[0] = endCenter[0]; onPlane[1] = endCenter[1]; onPlane[2] = endCenter[2]; if(! (lastOnPositive ^ secondLastOnPositive) ) { /* endId and endId+1 lie on the same side of the plane, so we need need endId-1 to calculate the intersection with the planar figure. this should exist since we know that the fiber crosses the planar figure*/ p = fiberPolyData->GetPoint( pointsInCell[ endId-1 ] ); p1[0] = p[0]; p1[1] = p[1]; p1[2] = p[2]; /* endId and endId+1 lie on the same side of the plane so startId is also part of the ROI*/ double *end = fiberPolyData->GetPoint( pointsInCell[endId] ); points->InsertNextPoint(end); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } d = ( (onPlane-p0)*normal) / ( (p0-p1) * normal ); newPoint = (p0-p1); newPoint[0] = d*newPoint[0] + p0[0]; newPoint[1] = d*newPoint[1] + p0[1]; newPoint[2] = d*newPoint[2] + p0[2]; insertPoint[0] = newPoint[0]; insertPoint[1] = newPoint[1]; insertPoint[2] = newPoint[2]; //Insert the Last Point (intersection with the end roi) points->InsertNextPoint(insertPoint); polyLine->GetPointIds()->InsertId(lineIndex,pointIndex); lineIndex++; pointIndex++; } // add polyline to vtk cell array cells->InsertNextCell(polyLine); } // Add the points to the dataset polyData->SetPoints(points); // Add the lines to the dataset polyData->SetLines(cells); mitk::FiberBundleX::Pointer cutBundle = mitk::FiberBundleX::New(polyData); mitk::DataNode::Pointer cutNode = mitk::DataNode::New(); cutNode->SetData(cutBundle); std::string name = "fiberCut"; cutNode->SetName(name); GetDataStorage()->Add(cutNode); } void QmitkTractbasedSpatialStatisticsView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkTractbasedSpatialStatisticsView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } -/* -void QmitkTractbasedSpatialStatisticsView::AdjustPlotMeasure(const QString & text) -{ - berry::ISelection::ConstPointer sel( - this->GetSite()->GetWorkbenchWindow()->GetSelectionService()->GetSelection("org.mitk.views.datamanager")); - m_CurrentSelection = sel.Cast(); - m_SelListener.Cast()->DoSelectionChanged(sel); -}*/ - - void QmitkTractbasedSpatialStatisticsView::CreateRoi() { - - // It is important to load the MeanFASkeletonMask image in MITK to make sure that point selection and - // pathfinding is done on the same image - //string filename = m_TbssWorkspaceManager.GetInputDir().toStdString() + "/stats/" + m_TbssWorkspaceManager.GetMeanFASkeletonMask().toStdString(); - - // Implement a way to obtain skeleton and skeletonFA without sml workspace - bool ok; double threshold = QInputDialog::getDouble(m_Controls->m_CreateRoi, tr("Set an FA threshold"), tr("Threshold:"), 0.2, 0.0, 1.0, 2, &ok); if(!ok) return; mitk::Image::Pointer image; - std::vector nodes = this->GetDataManagerSelection(); for ( int i=0; iGetData()->GetNameOfClass())==0) { mitk::Image* img = static_cast(nodes[i]->GetData()); if(img->GetDimension() == 3) { image = img; } } } - - if(image.IsNull()) { return; } mitk::TractAnalyzer analyzer; analyzer.SetInputImage(image); analyzer.SetThreshold(threshold); // Set Pointset to analyzer analyzer.SetPointSet(m_PointSetNode); // Run Analyzer - analyzer.BuildGraph(); + analyzer.MakeRoi(); // Obtain tbss roi image from analyzer mitk::TbssRoiImage::Pointer tbssRoi = analyzer.GetRoiImage(); tbssRoi->SetStructure(m_Controls->m_Structure->text().toStdString()); // get path description and set to interface std::string pathDescription = analyzer.GetPathDescription(); m_Controls->m_PathTextEdit->setPlainText(QString(pathDescription.c_str())); // Add roi image to datastorage AddTbssToDataStorage(tbssRoi, m_Controls->m_RoiName->text().toStdString()); } void QmitkTractbasedSpatialStatisticsView:: PlotFiberBundle(mitk::FiberBundleX *fib, mitk::Image* img, mitk::PlanarFigure* startRoi, mitk::PlanarFigure* endRoi) { bool avg = m_Controls->m_Average->isChecked(); int segments = m_Controls->m_Segments->value(); m_Controls->m_RoiPlotWidget->PlotFiberBetweenRois(fib, img, startRoi ,endRoi, avg, segments); m_Controls->m_RoiPlotWidget->SetPlottingFiber(true); mitk::RenderingManager::GetInstance()->ForceImmediateUpdateAll(); } void QmitkTractbasedSpatialStatisticsView::Plot(mitk::TbssImage* image, mitk::TbssRoiImage* roiImage) { if(m_Controls->m_TabWidget->currentWidget() == m_Controls->m_MeasureTAB) { std::vector< itk::Index<3> > roi = roiImage->GetRoi(); m_Roi = roi; m_CurrentGeometry = image->GetGeometry(); std::string resultfile = ""; std::string structure = roiImage->GetStructure(); m_Controls->m_RoiPlotWidget->SetGroups(image->GetGroupInfo()); m_Controls->m_RoiPlotWidget->SetProjections(image->GetImage()); m_Controls->m_RoiPlotWidget->SetRoi(roi); m_Controls->m_RoiPlotWidget->SetStructure(structure); m_Controls->m_RoiPlotWidget->SetMeasure( image->GetMeasurementInfo() ); m_Controls->m_RoiPlotWidget->DrawProfiles(resultfile); } m_Controls->m_RoiPlotWidget->SetPlottingFiber(false); } - - - - -void QmitkTractbasedSpatialStatisticsView::Masking() -{ - //QString filename = m_Controls->m_WorkingDirectory->text(); - QString filename = "E:/Experiments/tbss"; - QString faFiles = filename + "/AxD"; - QString maskFiles = filename + "/bin_masks"; - - - QDirIterator faDirIt(faFiles, QDir::Files | QDir::NoSymLinks, QDirIterator::Subdirectories); - QDirIterator maskDirIt(maskFiles, QDir::Files | QDir::NoSymLinks, QDirIterator::Subdirectories); - - std::vector faFilenames; - std::vector maskFilenames; - std::vector outputFilenames; - - while(faDirIt.hasNext() && maskDirIt.hasNext()) - { - faDirIt.next(); - maskDirIt.next(); - if((faDirIt.fileInfo().completeSuffix() == "nii" - || faDirIt.fileInfo().completeSuffix() == "mhd" - || faDirIt.fileInfo().completeSuffix() == "nii.gz") - && (maskDirIt.fileInfo().completeSuffix() == "nii" - || maskDirIt.fileInfo().completeSuffix() == "mhd" - || maskDirIt.fileInfo().completeSuffix() == "nii.gz")) - { - faFilenames.push_back(faDirIt.filePath().toStdString()); - outputFilenames.push_back(faDirIt.fileName().toStdString()); - maskFilenames.push_back(maskDirIt.filePath().toStdString()); - } - } - - std::vector::iterator faIt = faFilenames.begin(); - std::vector::iterator maskIt = maskFilenames.begin(); - std::vector::iterator outputIt = outputFilenames.begin(); - - // Now multiply all FA images with their corresponding masks - - QString outputDir = filename; - while(faIt != faFilenames.end() && maskIt != maskFilenames.end() && outputIt != outputFilenames.end()) - { - std::cout << "Mask " << *faIt << " with " << *maskIt << std::endl; - - typedef itk::MultiplyImageFilter MultiplicationFilterType; - - FloatReaderType::Pointer floatReader = FloatReaderType::New(); - CharReaderType::Pointer charReader = CharReaderType::New(); - - floatReader->SetFileName(*faIt); - //floatReader->Update(); - //FloatImageType::Pointer faImage = floatReader->GetOutput(); - - charReader->SetFileName(*maskIt); - //charReader->Update(); - // CharImageType::Pointer maskImage = charReader->GetOutput(); - - MultiplicationFilterType::Pointer multiplicationFilter = MultiplicationFilterType::New(); - multiplicationFilter->SetInput1(floatReader->GetOutput()); - multiplicationFilter->SetInput2(charReader->GetOutput()); - multiplicationFilter->Update(); - - //FloatImageType::Pointer maskedImage = FloatImageType::New(); - //maskedImage = MultiplicationFilter->GetOutput(); - - FloatWriterType::Pointer floatWriter = FloatWriterType::New(); - std::string s = faFiles.toStdString().append("/"+*outputIt); - floatWriter->SetFileName(s.c_str()); - floatWriter->SetInput(multiplicationFilter->GetOutput()); - floatWriter->Update(); - - ++faIt; - ++maskIt; - ++outputIt; - } -} 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 5ba0a0ebc7..0f3a45cfff 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkTractbasedSpatialStatisticsView.h @@ -1,255 +1,170 @@ /*=================================================================== 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 -#include - - #include #include "ui_QmitkTractbasedSpatialStatisticsViewControls.h" #include -#include -#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::Vector IntVectorType; -//typedef itk::VectorImage DirectionImageType; 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; -struct TbssSelListener; - - /*! - \brief QmitkTractbasedSpatialStatisticsView - - \warning This application module is not yet documented. Use "svn blame/praise/annotate" and ask the author to provide basic documentation. - - \sa QmitkFunctionalitymitkTbssWorkspaceManager - \ingroup Functionalities + * \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 { - - friend struct TbssSelListener; - - - - // this is needed for all Qt objesetupUicts that should have a Qt meta-object - // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkTractbasedSpatialStatisticsView(); virtual ~QmitkTractbasedSpatialStatisticsView(); virtual void CreateQtPartControl(QWidget *parent); /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); /// \brief Called when the functionality is activated virtual void Activated(); virtual void Deactivated(); protected slots: - - - void Masking(); - - - + // 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 InitPointsets(); + void PlotFiberBundle(mitk::FiberBundleX* fib, mitk::Image* img, mitk::PlanarFigure* startRoi=NULL, mitk::PlanarFigure* endRoi=NULL); - void SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name); - bool m_IsInitialized; + // 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; - - std::vector SortPoints(CharImageType::Pointer roi, CharImageType::IndexType currentPoint); - - bool PointVisited(std::vector points, CharImageType::IndexType point); - - // Modifies the current point by reference and returns true if no more points need to be visited - CharImageType::IndexType FindNextPoint(std::vector pointsVisited, - CharImageType::IndexType currentPoint, CharImageType::Pointer roi, bool &ready); - - - - - //void DoInitializeGridByVectorImage(FloatVectorImageType::Pointer vectorpic, CharImageType::Pointer roi ,std::string name); - - - - // Tokenizer needed for the roi files - void Tokenize(const std::string& str, - std::vector& tokens, - const std::string& delimiters = " ") - { - // Skip delimiters at beginning. - std::string::size_type lastPos = str.find_first_not_of(delimiters, 0); - // Find first "non-delimiter". - std::string::size_type pos = str.find_first_of(delimiters, lastPos); - - while (std::string::npos != pos || std::string::npos != lastPos) - { - // Found a token, add it to the vector. - tokens.push_back(str.substr(lastPos, pos - lastPos)); - // Skip delimiters. Note the "not_of" - lastPos = str.find_first_not_of(delimiters, pos); - // Find next "non-delimiter" - pos = str.find_first_of(delimiters, lastPos); - } - } - - mitk::DataNode::Pointer readNode(std::string f) - { - mitk::DataNode::Pointer node; - mitk::DataNodeFactory::Pointer nodeReader = mitk::DataNodeFactory::New(); - try - { - nodeReader->SetFileName(f); - nodeReader->Update(); - node = nodeReader->GetOutput(); - } - catch(...) { - MITK_ERROR << "Could not read file"; - return NULL; - } - - return node; - } - - /*template < typename TPixel, unsigned int VImageDimension > - void ToITK4D( itk::Image* inputImage, Float4DImageType::Pointer& outputImage );*/ - - - std::string ReadFile(std::string whatfile); - + // Used to save the region of interest in a vector of itk::index. std::vector< itk::Index<3> > m_Roi; mitk::FiberBundleX* m_Fib; - std::string m_CurrentStructure; - 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::TbssImage::Pointer m_CurrentTbssMetaImage; - - VectorImageType::Pointer ConvertToVectorImage(mitk::Image::Pointer mitkImg); - - mitk::DataNode::Pointer m_CurrentFiberNode; // needed for the index property when interacting with the plot widget - mitk::DataNode::Pointer m_CurrentStartRoi; // needed when a plot should only show values between a start end end roi - mitk::DataNode::Pointer m_CurrentEndRoi; // idem dito - + // 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