diff --git a/Modules/DiffusionImaging/Algorithms/itkMergeDiffusionImagesFilter.h b/Modules/DiffusionImaging/Algorithms/itkMergeDiffusionImagesFilter.h new file mode 100644 index 0000000000..e1706e5d5b --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkMergeDiffusionImagesFilter.h @@ -0,0 +1,100 @@ +/*=================================================================== + +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. + +===================================================================*/ +/*========================================================================= + +Program: Tensor ToolKit - TTK +Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkMergeDiffusionImagesFilter.h $ +Language: C++ +Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ +Version: $Revision: 68 $ + +Copyright (c) INRIA 2010. All rights reserved. +See LICENSE.txt for details. + +This software is distributed WITHOUT ANY WARRANTY; without even +the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR +PURPOSE. See the above copyright notices for more information. + +=========================================================================*/ +#ifndef _itk_MergeDiffusionImagesFilter_h_ +#define _itk_MergeDiffusionImagesFilter_h_ + +#include "itkImageToImageFilter.h" +#include +#include +#include + +namespace itk +{ + +template +class MergeDiffusionImagesFilter + : public ImageSource > +{ + +public: + + typedef itk::VectorImage DwiImageType; + typedef MergeDiffusionImagesFilter Self; + typedef ImageSource Superclass; + + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + typedef typename DwiImageType::PixelType DwiPixelType; + typedef typename DwiImageType::RegionType DwiRegionType; + typedef typename std::vector< typename DwiImageType::Pointer > DwiImageContainerType; + + typedef vnl_vector_fixed< double, 3 > GradientType; + typedef itk::VectorContainer< unsigned int, GradientType > GradientListType; + typedef typename std::vector< GradientListType::Pointer > GradientListContainerType; + + itkTypeMacro (MergeDiffusionImagesFilter, ImageToImageFilter) + itkStaticConstMacro (ImageDimension, unsigned int, DwiImageType::ImageDimension); + itkNewMacro (Self) + + void SetImageVolumes(DwiImageContainerType cont){ m_ImageVolumes=cont; } + void SetGradientLists(GradientListContainerType cont){ m_GradientLists=cont; } + void SetBValues(std::vector< double > cont){ m_BValues=cont; } + + itkGetMacro( OutputGradients, GradientListType::Pointer ) + itkGetMacro( BValue, double ) + +protected: + + MergeDiffusionImagesFilter(){} + ~MergeDiffusionImagesFilter(){} + + void GenerateData(); + + DwiImageContainerType m_ImageVolumes; + GradientListContainerType m_GradientLists; + std::vector< double > m_BValues; + int m_NumGradients; + GradientListType::Pointer m_OutputGradients; + double m_BValue; +}; + + +} // end of namespace + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkMergeDiffusionImagesFilter.txx" +#endif + + +#endif diff --git a/Modules/DiffusionImaging/Algorithms/itkMergeDiffusionImagesFilter.txx b/Modules/DiffusionImaging/Algorithms/itkMergeDiffusionImagesFilter.txx new file mode 100644 index 0000000000..9dd0e145f1 --- /dev/null +++ b/Modules/DiffusionImaging/Algorithms/itkMergeDiffusionImagesFilter.txx @@ -0,0 +1,128 @@ +/*=================================================================== + +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. + +===================================================================*/ +/*========================================================================= + 2 + 3 Program: Tensor ToolKit - TTK + 4 Module: $URL: svn://scm.gforge.inria.fr/svn/ttk/trunk/Algorithms/itkMergeDiffusionImagesFilter.txx $ + 5 Language: C++ + 6 Date: $Date: 2010-06-07 13:39:13 +0200 (Mo, 07 Jun 2010) $ + 7 Version: $Revision: 68 $ + 8 + 9 Copyright (c) INRIA 2010. All rights reserved. + 10 See LICENSE.txt for details. + 11 + 12 This software is distributed WITHOUT ANY WARRANTY; without even + 13 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + 14 PURPOSE. See the above copyright notices for more information. + 15 + 16 =========================================================================*/ +#ifndef _itk_MergeDiffusionImagesFilter_txx_ +#define _itk_MergeDiffusionImagesFilter_txx_ +#endif + +#include "itkMergeDiffusionImagesFilter.h" +#include "itkTensorToL2NormImageFilter.h" +#include "itkRescaleIntensityImageFilter.h" +#include +#include +#include + +namespace itk +{ + +template +void +MergeDiffusionImagesFilter +::GenerateData () +{ + + if( m_ImageVolumes.size()<2 ) + throw itk::ExceptionObject (__FILE__,__LINE__,"Error: cannot combine less than two DWIs."); + + if( m_GradientLists.size()!=m_ImageVolumes.size() || m_ImageVolumes.size()!=m_BValues.size() || m_BValues.size()!=m_GradientLists.size() ) + throw itk::ExceptionObject (__FILE__,__LINE__,"Error: need same number of b-values, image volumes and gradient containers."); + + typename DwiImageType::Pointer img = m_ImageVolumes.at(0); + + m_NumGradients = 0; + for (int i=0; iSize(); + typename DwiImageType::Pointer tmp = m_ImageVolumes.at(i); + if ( img->GetLargestPossibleRegion()!=tmp->GetLargestPossibleRegion() ) + throw itk::ExceptionObject (__FILE__,__LINE__,"Error: images are not of same size."); + } + + m_BValue = m_BValues.at(0); + m_OutputGradients = GradientListType::New(); + + typename DwiImageType::Pointer outImage = DwiImageType::New(); + outImage->SetSpacing( img->GetSpacing() ); // Set the image spacing + outImage->SetOrigin( img->GetOrigin() ); // Set the image origin + outImage->SetDirection( img->GetDirection() ); // Set the image direction + outImage->SetLargestPossibleRegion( img->GetLargestPossibleRegion()); + outImage->SetBufferedRegion( img->GetLargestPossibleRegion() ); + outImage->SetRequestedRegion( img->GetLargestPossibleRegion() ); + outImage->SetVectorLength(m_NumGradients); + outImage->Allocate(); + + this->SetNumberOfRequiredOutputs(1); + this->SetNthOutput (0, outImage); + + typedef ImageRegionIterator IteratorOutputType; + IteratorOutputType itOut (this->GetOutput(0), this->GetOutput(0)->GetLargestPossibleRegion()); + + MITK_INFO << "MergeDiffusionImagesFilter: merging images"; + GradientType zeroG; zeroG.fill(0.0); + boost::progress_display disp(this->GetOutput(0)->GetLargestPossibleRegion().GetNumberOfPixels()); + while(!itOut.IsAtEnd()) + { + ++disp; + DwiPixelType out; + out.SetSize(m_NumGradients); + out.Fill(0); + + int c=0; + for (int i=0; iSize(); j++) + { + GradientType g = gradients->GetElement(j); + + if (g.magnitude()>0.0001) + { + g.normalize(); + g *= std::sqrt(bValue/m_BValue); + } + else + g = zeroG; + + m_OutputGradients->InsertElement(c, g); + out[c] = static_cast(img->GetPixel(itOut.GetIndex())[j]); + c++; + } + } + + itOut.Set(out); + ++itOut; + } +} + +} // end of namespace diff --git a/Modules/DiffusionImaging/Algorithms/itkTractDensityImageFilter.cpp b/Modules/DiffusionImaging/Algorithms/itkTractDensityImageFilter.cpp index 1902bb1613..39d27bcdd2 100644 --- a/Modules/DiffusionImaging/Algorithms/itkTractDensityImageFilter.cpp +++ b/Modules/DiffusionImaging/Algorithms/itkTractDensityImageFilter.cpp @@ -1,223 +1,226 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Coindex[1]right (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 "itkTractDensityImageFilter.h" // VTK #include #include #include // misc #include +#include namespace itk{ template< class OutputImageType > TractDensityImageFilter< OutputImageType >::TractDensityImageFilter() : m_BinaryOutput(false) , m_InvertImage(false) , m_UpsamplingFactor(1) , m_InputImage(NULL) , m_UseImageGeometry(false) , m_OutputAbsoluteValues(false) { } template< class OutputImageType > TractDensityImageFilter< OutputImageType >::~TractDensityImageFilter() { } template< class OutputImageType > itk::Point TractDensityImageFilter< OutputImageType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class OutputImageType > void TractDensityImageFilter< OutputImageType >::GenerateData() { // generate upsampled image mitk::Geometry3D::Pointer geometry = m_FiberBundle->GetGeometry(); typename OutputImageType::Pointer outImage = this->GetOutput(); // calculate new image parameters mitk::Vector3D newSpacing; mitk::Point3D newOrigin; itk::Matrix newDirection; ImageRegion<3> upsampledRegion; if (m_UseImageGeometry && !m_InputImage.IsNull()) { MITK_INFO << "TractDensityImageFilter: using image geometry"; newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor; upsampledRegion = m_InputImage->GetLargestPossibleRegion(); newOrigin = m_InputImage->GetOrigin(); typename OutputImageType::RegionType::SizeType size = upsampledRegion.GetSize(); size[0] *= m_UpsamplingFactor; size[1] *= m_UpsamplingFactor; size[2] *= m_UpsamplingFactor; upsampledRegion.SetSize(size); newDirection = m_InputImage->GetDirection(); } else { MITK_INFO << "TractDensityImageFilter: using fiber bundle geometry"; newSpacing = geometry->GetSpacing()/m_UpsamplingFactor; newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); for (int i=0; i<3; i++) for (int j=0; j<3; j++) newDirection[j][i] = geometry->GetMatrixColumn(i)[j]; upsampledRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); upsampledRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); upsampledRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); } typename OutputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize(); // apply new image parameters outImage->SetSpacing( newSpacing ); outImage->SetOrigin( newOrigin ); outImage->SetDirection( newDirection ); outImage->SetRegions( upsampledRegion ); outImage->Allocate(); outImage->FillBuffer(0.0); int w = upsampledSize[0]; int h = upsampledSize[1]; int d = upsampledSize[2]; // set/initialize output OutPixelType* outImageBufferPointer = (OutPixelType*)outImage->GetBufferPointer(); // resample fiber bundle float minSpacing = 1; if(newSpacing[0]GetDeepCopy(); m_FiberBundle->ResampleFibers(minSpacing); MITK_INFO << "TractDensityImageFilter: starting image generation"; vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int numFibers = m_FiberBundle->GetNumFibers(); + boost::progress_display disp(numFibers); for( int i=0; iGetNextCell ( numPoints, points ); // fill output image for( int j=0; j vertex = GetItkPoint(fiberPolyData->GetPoint(points[j])); itk::Index<3> index; itk::ContinuousIndex contIndex; outImage->TransformPhysicalPointToIndex(vertex, index); outImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex); float frac_x = contIndex[0] - index[0]; float frac_y = contIndex[1] - index[1]; float frac_z = contIndex[2] - index[2]; if (frac_x<0) { index[0] -= 1; frac_x += 1; } if (frac_y<0) { index[1] -= 1; frac_y += 1; } if (frac_z<0) { index[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (index[0] < 0 || index[0] >= w-1) continue; if (index[1] < 0 || index[1] >= h-1) continue; if (index[2] < 0 || index[2] >= d-1) continue; if (m_BinaryOutput) { outImageBufferPointer[( index[0] + w*(index[1] + h*index[2] ))] = 1; outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2] ))] = 1; outImageBufferPointer[( index[0] + w*(index[1] + h*index[2]+h))] = 1; outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2]+h))] = 1; outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2] ))] = 1; outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2]+h))] = 1; outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2] ))] = 1; outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2]+h))] = 1; } else { outImageBufferPointer[( index[0] + w*(index[1] + h*index[2] ))] += ( frac_x)*( frac_y)*( frac_z); outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2] ))] += ( frac_x)*(1-frac_y)*( frac_z); outImageBufferPointer[( index[0] + w*(index[1] + h*index[2]+h))] += ( frac_x)*( frac_y)*(1-frac_z); outImageBufferPointer[( index[0] + w*(index[1]+1+ h*index[2]+h))] += ( frac_x)*(1-frac_y)*(1-frac_z); outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2] ))] += (1-frac_x)*( frac_y)*( frac_z); outImageBufferPointer[( index[0]+1 + w*(index[1] + h*index[2]+h))] += (1-frac_x)*( frac_y)*(1-frac_z); outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2] ))] += (1-frac_x)*(1-frac_y)*( frac_z); outImageBufferPointer[( index[0]+1 + w*(index[1]+1+ h*index[2]+h))] += (1-frac_x)*(1-frac_y)*(1-frac_z); } } } if (!m_OutputAbsoluteValues && !m_BinaryOutput) { MITK_INFO << "TractDensityImageFilter: max-normalizing output image"; OutPixelType max = 0; for (int i=0; i0) for (int i=0; i #include #include +#include namespace itk{ template< class OutputImageType > TractsToFiberEndingsImageFilter< OutputImageType >::TractsToFiberEndingsImageFilter() : m_InvertImage(false) , m_UpsamplingFactor(1) , m_InputImage(NULL) , m_UseImageGeometry(false) { } template< class OutputImageType > TractsToFiberEndingsImageFilter< OutputImageType >::~TractsToFiberEndingsImageFilter() { } template< class OutputImageType > itk::Point TractsToFiberEndingsImageFilter< OutputImageType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class OutputImageType > void TractsToFiberEndingsImageFilter< OutputImageType >::GenerateData() { // generate upsampled image mitk::Geometry3D::Pointer geometry = m_FiberBundle->GetGeometry(); typename OutputImageType::Pointer outImage = this->GetOutput(); // calculate new image parameters mitk::Vector3D newSpacing; mitk::Point3D newOrigin; itk::Matrix newDirection; ImageRegion<3> upsampledRegion; if (m_UseImageGeometry && !m_InputImage.IsNull()) { newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor; upsampledRegion = m_InputImage->GetLargestPossibleRegion(); newOrigin = m_InputImage->GetOrigin(); typename OutputImageType::RegionType::SizeType size = upsampledRegion.GetSize(); size[0] *= m_UpsamplingFactor; size[1] *= m_UpsamplingFactor; size[2] *= m_UpsamplingFactor; upsampledRegion.SetSize(size); newDirection = m_InputImage->GetDirection(); } else { newSpacing = geometry->GetSpacing()/m_UpsamplingFactor; newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); for (int i=0; i<3; i++) for (int j=0; j<3; j++) newDirection[j][i] = geometry->GetMatrixColumn(i)[j]; upsampledRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); upsampledRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); upsampledRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); } typename OutputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize(); // apply new image parameters outImage->SetSpacing( newSpacing ); outImage->SetOrigin( newOrigin ); outImage->SetDirection( newDirection ); outImage->SetRegions( upsampledRegion ); outImage->Allocate(); int w = upsampledSize[0]; int h = upsampledSize[1]; int d = upsampledSize[2]; // set/initialize output OutPixelType* outImageBufferPointer = (OutPixelType*)outImage->GetBufferPointer(); for (int i=0; i fiberPolyData = m_FiberBundle->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int numFibers = m_FiberBundle->GetNumFibers(); + boost::progress_display disp(numFibers); for( int i=0; iGetNextCell ( numPoints, points ); // fill output image if (numPoints>0) { itk::Point vertex = GetItkPoint(fiberPolyData->GetPoint(points[0])); itk::Index<3> index; outImage->TransformPhysicalPointToIndex(vertex, index); outImage->SetPixel(index, 1); } if (numPoints>2) { itk::Point vertex = GetItkPoint(fiberPolyData->GetPoint(points[numPoints-1])); itk::Index<3> index; outImage->TransformPhysicalPointToIndex(vertex, index); outImage->SetPixel(index, 1); } } if (m_InvertImage) for (int i=0; i #include #include // misc #include +#include namespace itk{ template< class OutputImageType > TractsToRgbaImageFilter< OutputImageType >::TractsToRgbaImageFilter() : m_InvertImage(false) , m_UpsamplingFactor(1) , m_InputImage(NULL) , m_UseImageGeometry(false) { } template< class OutputImageType > TractsToRgbaImageFilter< OutputImageType >::~TractsToRgbaImageFilter() { } template< class OutputImageType > itk::Point TractsToRgbaImageFilter< OutputImageType >::GetItkPoint(double point[3]) { itk::Point itkPoint; itkPoint[0] = point[0]; itkPoint[1] = point[1]; itkPoint[2] = point[2]; return itkPoint; } template< class OutputImageType > void TractsToRgbaImageFilter< OutputImageType >::GenerateData() { if(&typeid(OutPixelType) != &typeid(itk::RGBAPixel)) return; // generate upsampled image mitk::Geometry3D::Pointer geometry = m_FiberBundle->GetGeometry(); typename OutputImageType::Pointer outImage = this->GetOutput(); // calculate new image parameters mitk::Vector3D newSpacing; mitk::Point3D newOrigin; itk::Matrix newDirection; ImageRegion<3> upsampledRegion; if (m_UseImageGeometry && !m_InputImage.IsNull()) { newSpacing = m_InputImage->GetSpacing()/m_UpsamplingFactor; upsampledRegion = m_InputImage->GetLargestPossibleRegion(); newOrigin = m_InputImage->GetOrigin(); typename OutputImageType::RegionType::SizeType size = upsampledRegion.GetSize(); size[0] *= m_UpsamplingFactor; size[1] *= m_UpsamplingFactor; size[2] *= m_UpsamplingFactor; upsampledRegion.SetSize(size); newDirection = m_InputImage->GetDirection(); } else { newSpacing = geometry->GetSpacing()/m_UpsamplingFactor; newOrigin = geometry->GetOrigin(); mitk::Geometry3D::BoundsArrayType bounds = geometry->GetBounds(); newOrigin[0] += bounds.GetElement(0); newOrigin[1] += bounds.GetElement(2); newOrigin[2] += bounds.GetElement(4); for (int i=0; i<3; i++) for (int j=0; j<3; j++) newDirection[j][i] = geometry->GetMatrixColumn(i)[j]; upsampledRegion.SetSize(0, geometry->GetExtent(0)*m_UpsamplingFactor); upsampledRegion.SetSize(1, geometry->GetExtent(1)*m_UpsamplingFactor); upsampledRegion.SetSize(2, geometry->GetExtent(2)*m_UpsamplingFactor); } typename OutputImageType::RegionType::SizeType upsampledSize = upsampledRegion.GetSize(); // apply new image parameters outImage->SetSpacing( newSpacing ); outImage->SetOrigin( newOrigin ); outImage->SetDirection( newDirection ); outImage->SetRegions( upsampledRegion ); outImage->Allocate(); int w = upsampledSize[0]; int h = upsampledSize[1]; int d = upsampledSize[2]; // set/initialize output unsigned char* outImageBufferPointer = (unsigned char*)outImage->GetBufferPointer(); float* buffer = new float[w*h*d*4]; for (int i=0; iGetDeepCopy(); m_FiberBundle->ResampleFibers(minSpacing); vtkSmartPointer fiberPolyData = m_FiberBundle->GetFiberPolyData(); vtkSmartPointer vLines = fiberPolyData->GetLines(); vLines->InitTraversal(); int numFibers = m_FiberBundle->GetNumFibers(); + boost::progress_display disp(numFibers); for( int i=0; iGetNextCell ( numPoints, points ); // calc directions (which are used as weights) std::list< itk::Point > rgbweights; std::list intensities; for( int j=0; j vertex = GetItkPoint(fiberPolyData->GetPoint(points[j])); itk::Point vertexPost = GetItkPoint(fiberPolyData->GetPoint(points[j+1])); itk::Point dir; dir[0] = fabs((vertexPost[0] - vertex[0]) * outImage->GetSpacing()[0]); dir[1] = fabs((vertexPost[1] - vertex[1]) * outImage->GetSpacing()[1]); dir[2] = fabs((vertexPost[2] - vertex[2]) * outImage->GetSpacing()[2]); rgbweights.push_back(dir); float intensity = sqrt(dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2]); intensities.push_back(intensity); // last point gets same as previous one if(j==numPoints-2) { rgbweights.push_back(dir); intensities.push_back(intensity); } } // fill output image for( int j=0; j vertex = GetItkPoint(fiberPolyData->GetPoint(points[j])); itk::Index<3> index; itk::ContinuousIndex contIndex; outImage->TransformPhysicalPointToIndex(vertex, index); outImage->TransformPhysicalPointToContinuousIndex(vertex, contIndex); float frac_x = contIndex[0] - index[0]; float frac_y = contIndex[1] - index[1]; float frac_z = contIndex[2] - index[2]; int px = index[0]; if (frac_x<0) { px -= 1; frac_x += 1; } int py = index[1]; if (frac_y<0) { py -= 1; frac_y += 1; } int pz = index[2]; if (frac_z<0) { pz -= 1; frac_z += 1; } // int coordinates inside image? if (px < 0 || px >= w-1) continue; if (py < 0 || py >= h-1) continue; if (pz < 0 || pz >= d-1) continue; float scale = 100 * pow((float)m_UpsamplingFactor,3); itk::Point rgbweight = rgbweights.front(); rgbweights.pop_front(); float intweight = intensities.front(); intensities.pop_front(); // add to r-channel in output image buffer[0+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[0] * scale; buffer[0+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * rgbweight[0] * scale; buffer[0+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * rgbweight[0] * scale; buffer[0+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * rgbweight[0] * scale; buffer[0+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[0] * scale; buffer[0+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * rgbweight[0] * scale; buffer[0+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * rgbweight[0] * scale; buffer[0+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * rgbweight[0] * scale; // add to g-channel in output image buffer[1+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[1] * scale; buffer[1+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * rgbweight[1] * scale; buffer[1+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * rgbweight[1] * scale; buffer[1+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * rgbweight[1] * scale; buffer[1+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[1] * scale; buffer[1+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * rgbweight[1] * scale; buffer[1+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * rgbweight[1] * scale; buffer[1+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * rgbweight[1] * scale; // add to b-channel in output image buffer[2+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[2] * scale; buffer[2+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * rgbweight[2] * scale; buffer[2+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * rgbweight[2] * scale; buffer[2+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * rgbweight[2] * scale; buffer[2+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * rgbweight[2] * scale; buffer[2+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * rgbweight[2] * scale; buffer[2+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * rgbweight[2] * scale; buffer[2+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * rgbweight[2] * scale; // add to a-channel in output image buffer[3+4*( px + w*(py + h*pz ))] += (1-frac_x)*(1-frac_y)*(1-frac_z) * intweight * scale; buffer[3+4*( px + w*(py+1+ h*pz ))] += (1-frac_x)*( frac_y)*(1-frac_z) * intweight * scale; buffer[3+4*( px + w*(py + h*pz+h))] += (1-frac_x)*(1-frac_y)*( frac_z) * intweight * scale; buffer[3+4*( px + w*(py+1+ h*pz+h))] += (1-frac_x)*( frac_y)*( frac_z) * intweight * scale; buffer[3+4*( px+1 + w*(py + h*pz ))] += ( frac_x)*(1-frac_y)*(1-frac_z) * intweight * scale; buffer[3+4*( px+1 + w*(py + h*pz+h))] += ( frac_x)*(1-frac_y)*( frac_z) * intweight * scale; buffer[3+4*( px+1 + w*(py+1+ h*pz ))] += ( frac_x)*( frac_y)*(1-frac_z) * intweight * scale; buffer[3+4*( px+1 + w*(py+1+ h*pz+h))] += ( frac_x)*( frac_y)*( frac_z) * intweight * scale; } } float maxRgb = 0.000000001; float maxInt = 0.000000001; int numPix; numPix = w*h*d*4; // calc maxima for(int i=0; i maxRgb) maxRgb = buffer[i]; } else { if(buffer[i] > maxInt) maxInt = buffer[i]; } } // write output, normalized uchar 0..255 for(int i=0; i #include #include #include "itkStreamlineTrackingFilter.h" #include #include #include #define _USE_MATH_DEFINES #include namespace itk { //#define QBALL_RECON_PI M_PI template< class TTensorPixelType, class TPDPixelType> StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::StreamlineTrackingFilter(): m_FaThreshold(0.2), m_StepSize(1), m_MaxLength(10000), m_SeedsPerVoxel(1), m_F(1.0), m_G(0.0), m_Interpolate(true), m_MinTractLength(0.0) { // At least 1 inputs is necessary for a vector image. // For images added one at a time we need at least six this->SetNumberOfRequiredInputs( 1 ); } template< class TTensorPixelType, class TPDPixelType> double StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::RoundToNearest(double num) { return (num > 0.0) ? floor(num + 0.5) : ceil(num - 0.5); } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::BeforeThreadedGenerateData() { m_FiberPolyData = FiberPolyDataType::New(); m_Points = vtkPoints::New(); m_Cells = vtkCellArray::New(); m_InputImage = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); m_ImageSize.resize(3); m_ImageSize[0] = m_InputImage->GetLargestPossibleRegion().GetSize()[0]; m_ImageSize[1] = m_InputImage->GetLargestPossibleRegion().GetSize()[1]; m_ImageSize[2] = m_InputImage->GetLargestPossibleRegion().GetSize()[2]; if (m_ImageSize[0]<3 || m_ImageSize[1]<3 || m_ImageSize[2]<3) m_Interpolate = false; m_ImageSpacing.resize(3); m_ImageSpacing[0] = m_InputImage->GetSpacing()[0]; m_ImageSpacing[1] = m_InputImage->GetSpacing()[1]; m_ImageSpacing[2] = m_InputImage->GetSpacing()[2]; float minSpacing; if(m_ImageSpacing[0]::New(); for (int i=0; iGetNumberOfThreads(); i++) { FiberPolyDataType poly = FiberPolyDataType::New(); m_PolyDataContainer->InsertElement(i, poly); } if (m_SeedImage.IsNull()) { // initialize mask image m_SeedImage = ItkUcharImgType::New(); m_SeedImage->SetSpacing( m_InputImage->GetSpacing() ); m_SeedImage->SetOrigin( m_InputImage->GetOrigin() ); m_SeedImage->SetDirection( m_InputImage->GetDirection() ); m_SeedImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_SeedImage->Allocate(); m_SeedImage->FillBuffer(1); } if (m_MaskImage.IsNull()) { // initialize mask image m_MaskImage = ItkUcharImgType::New(); m_MaskImage->SetSpacing( m_InputImage->GetSpacing() ); m_MaskImage->SetOrigin( m_InputImage->GetOrigin() ); m_MaskImage->SetDirection( m_InputImage->GetDirection() ); m_MaskImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1); } m_FaImage = ItkFloatImgType::New(); m_FaImage->SetSpacing( m_InputImage->GetSpacing() ); m_FaImage->SetOrigin( m_InputImage->GetOrigin() ); m_FaImage->SetDirection( m_InputImage->GetDirection() ); m_FaImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_FaImage->Allocate(); m_FaImage->FillBuffer(0.0); m_PdImage = ItkPDImgType::New(); m_PdImage->SetSpacing( m_InputImage->GetSpacing() ); m_PdImage->SetOrigin( m_InputImage->GetOrigin() ); m_PdImage->SetDirection( m_InputImage->GetDirection() ); m_PdImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_PdImage->Allocate(); m_EmaxImage = ItkFloatImgType::New(); m_EmaxImage->SetSpacing( m_InputImage->GetSpacing() ); m_EmaxImage->SetOrigin( m_InputImage->GetOrigin() ); m_EmaxImage->SetDirection( m_InputImage->GetDirection() ); m_EmaxImage->SetRegions( m_InputImage->GetLargestPossibleRegion() ); m_EmaxImage->Allocate(); m_EmaxImage->FillBuffer(1.0); typedef itk::DiffusionTensor3D TensorType; typename TensorType::EigenValuesArrayType eigenvalues; typename TensorType::EigenVectorsMatrixType eigenvectors; for (int x=0; xGetPixel(index); vnl_vector_fixed dir; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); dir.normalize(); m_PdImage->SetPixel(index, dir); m_FaImage->SetPixel(index, tensor.GetFractionalAnisotropy()); m_EmaxImage->SetPixel(index, 2/eigenvalues[2]); } if (m_Interpolate) std::cout << "StreamlineTrackingFilter: using trilinear interpolation" << std::endl; else { if (m_MinCurvatureRadius<0.0) m_MinCurvatureRadius = 0.1*minSpacing; std::cout << "StreamlineTrackingFilter: using nearest neighbor interpolation" << std::endl; } if (m_MinCurvatureRadius<0.0) m_MinCurvatureRadius = 0.5*minSpacing; std::cout << "StreamlineTrackingFilter: Min. curvature radius: " << m_MinCurvatureRadius << std::endl; std::cout << "StreamlineTrackingFilter: FA threshold: " << m_FaThreshold << std::endl; std::cout << "StreamlineTrackingFilter: stepsize: " << m_StepSize << " mm" << std::endl; std::cout << "StreamlineTrackingFilter: f: " << m_F << std::endl; std::cout << "StreamlineTrackingFilter: g: " << m_G << std::endl; std::cout << "StreamlineTrackingFilter: starting streamline tracking" << std::endl; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::CalculateNewPosition(itk::ContinuousIndex& pos, vnl_vector_fixed& dir, typename InputImageType::IndexType& index) { if (true) { dir *= m_StepSize; pos[0] += dir[0]/m_ImageSpacing[0]; pos[1] += dir[1]/m_ImageSpacing[1]; pos[2] += dir[2]/m_ImageSpacing[2]; index[0] = RoundToNearest(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); } else { dir[0] /= m_ImageSpacing[0]; dir[1] /= m_ImageSpacing[1]; dir[2] /= m_ImageSpacing[2]; int smallest = 0; float x = 100000; if (dir[0]>0) { if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps) x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0]; else x = fabs(pos[0]-std::ceil(pos[0])-0.5)/dir[0]; } else if (dir[0]<0) { if (fabs(fabs(RoundToNearest(pos[0])-pos[0])-0.5)>mitk::eps) x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0]; else x = -fabs(pos[0]-std::floor(pos[0])+0.5)/dir[0]; } float s = x; float y = 100000; if (dir[1]>0) { if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps) y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1]; else y = fabs(pos[1]-std::ceil(pos[1])-0.5)/dir[1]; } else if (dir[1]<0) { if (fabs(fabs(RoundToNearest(pos[1])-pos[1])-0.5)>mitk::eps) y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1]; else y = -fabs(pos[1]-std::floor(pos[1])+0.5)/dir[1]; } if (s>y) { s=y; smallest = 1; } float z = 100000; if (dir[2]>0) { if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps) z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2]; else z = fabs(pos[2]-std::ceil(pos[2])-0.5)/dir[2]; } else if (dir[2]<0) { if (fabs(fabs(RoundToNearest(pos[2])-pos[2])-0.5)>mitk::eps) z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2]; else z = -fabs(pos[2]-std::floor(pos[2])+0.5)/dir[2]; } if (s>z) { s=z; smallest = 2; } // MITK_INFO << "---------------------------------------------"; // MITK_INFO << "s: " << s; // MITK_INFO << "dir: " << dir; // MITK_INFO << "old: " << pos[0] << ", " << pos[1] << ", " << pos[2]; pos[0] += dir[0]*s; pos[1] += dir[1]*s; pos[2] += dir[2]*s; switch (smallest) { case 0: if (dir[0]<0) index[0] = std::floor(pos[0]); else index[0] = std::ceil(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); break; case 1: if (dir[1]<0) index[1] = std::floor(pos[1]); else index[1] = std::ceil(pos[1]); index[0] = RoundToNearest(pos[0]); index[2] = RoundToNearest(pos[2]); break; case 2: if (dir[2]<0) index[2] = std::floor(pos[2]); else index[2] = std::ceil(pos[2]); index[1] = RoundToNearest(pos[1]); index[0] = RoundToNearest(pos[0]); } // float x = 100000; // if (dir[0]>0) // x = fabs(pos[0]-RoundToNearest(pos[0])-0.5)/dir[0]; // else if (dir[0]<0) // x = -fabs(pos[0]-RoundToNearest(pos[0])+0.5)/dir[0]; // float s = x; // float y = 100000; // if (dir[1]>0) // y = fabs(pos[1]-RoundToNearest(pos[1])-0.5)/dir[1]; // else if (dir[1]<0) // y = -fabs(pos[1]-RoundToNearest(pos[1])+0.5)/dir[1]; // if (s>y) // s=y; // float z = 100000; // if (dir[2]>0) // z = fabs(pos[2]-RoundToNearest(pos[2])-0.5)/dir[2]; // else if (dir[2]<0) // z = -fabs(pos[2]-RoundToNearest(pos[2])+0.5)/dir[2]; // if (s>z) // s=z; // s *= 1.001; // pos[0] += dir[0]*s; // pos[1] += dir[1]*s; // pos[2] += dir[2]*s; // index[0] = RoundToNearest(pos[0]); // index[1] = RoundToNearest(pos[1]); // index[2] = RoundToNearest(pos[2]); // MITK_INFO << "new: " << pos[0] << ", " << pos[1] << ", " << pos[2]; } } template< class TTensorPixelType, class TPDPixelType> bool StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::IsValidPosition(itk::ContinuousIndex& pos, typename InputImageType::IndexType &index, vnl_vector_fixed< float, 8 >& interpWeights) { - if (m_MaskImage->GetPixel(index)==0 || !m_InputImage->GetLargestPossibleRegion().IsInside(index)) + if (!m_InputImage->GetLargestPossibleRegion().IsInside(index) || m_MaskImage->GetPixel(index)==0) return false; if (m_Interpolate) { float frac_x = pos[0] - index[0]; float frac_y = pos[1] - index[1]; float frac_z = pos[2] - index[2]; if (frac_x<0) { index[0] -= 1; frac_x += 1; } if (frac_y<0) { index[1] -= 1; frac_y += 1; } if (frac_z<0) { index[2] -= 1; frac_z += 1; } frac_x = 1-frac_x; frac_y = 1-frac_y; frac_z = 1-frac_z; // int coordinates inside image? if (index[0] < 0 || index[0] >= m_ImageSize[0]-1) return false; if (index[1] < 0 || index[1] >= m_ImageSize[1]-1) return false; if (index[2] < 0 || index[2] >= m_ImageSize[2]-1) return false; interpWeights[0] = ( frac_x)*( frac_y)*( frac_z); interpWeights[1] = (1-frac_x)*( frac_y)*( frac_z); interpWeights[2] = ( frac_x)*(1-frac_y)*( frac_z); interpWeights[3] = ( frac_x)*( frac_y)*(1-frac_z); interpWeights[4] = (1-frac_x)*(1-frac_y)*( frac_z); interpWeights[5] = ( frac_x)*(1-frac_y)*(1-frac_z); interpWeights[6] = (1-frac_x)*( frac_y)*(1-frac_z); interpWeights[7] = (1-frac_x)*(1-frac_y)*(1-frac_z); typename InputImageType::IndexType tmpIdx; float FA = m_FaImage->GetPixel(index) * interpWeights[0]; tmpIdx = index; tmpIdx[0]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = index; tmpIdx[1]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = index; tmpIdx[2]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; FA += m_FaImage->GetPixel(tmpIdx) * interpWeights[7]; if (FAGetPixel(index) float StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::FollowStreamline(itk::ContinuousIndex pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids) { float tractLength = 0; typedef itk::DiffusionTensor3D TensorType; typename TensorType::EigenValuesArrayType eigenvalues; typename TensorType::EigenVectorsMatrixType eigenvectors; vnl_vector_fixed< float, 8 > interpWeights; typename InputImageType::IndexType index, indexOld; indexOld[0] = -1; indexOld[1] = -1; indexOld[2] = -1; itk::Point worldPos; float distance = 0; float distanceInVoxel = 0; // starting index and direction index[0] = RoundToNearest(pos[0]); index[1] = RoundToNearest(pos[1]); index[2] = RoundToNearest(pos[2]); vnl_vector_fixed dir = m_PdImage->GetPixel(index); dir *= dirSign; // reverse direction vnl_vector_fixed dirOld = dir; if (dir.magnitude()TransformContinuousIndexToPhysicalPoint( pos, worldPos ); ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); return tractLength; } else if (distance>=m_PointPistance) { m_InputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); distance = 0; } if (!m_Interpolate) // use nearest neighbour interpolation { if (indexOld!=index) // did we enter a new voxel? if yes, calculate new direction { dir = m_PdImage->GetPixel(index); // get principal direction typename InputImageType::PixelType tensor = m_InputImage->GetPixel(index); float scale = m_EmaxImage->GetPixel(index); dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2])); dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2])); dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2])); dir.normalize(); float angle = dot_product(dirOld, dir); if (angle<0) { dir *= -1; angle *= -1; } float r = m_StepSize/(2*std::asin(std::acos(angle)/2)); if (rGetPixel(index) * interpWeights[0]; typename InputImageType::IndexType tmpIdx = index; tmpIdx[0]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[1]; tmpIdx = index; tmpIdx[1]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = index; tmpIdx[2]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[7]; tensor.ComputeEigenAnalysis(eigenvalues, eigenvectors); dir[0] = eigenvectors(2, 0); dir[1] = eigenvectors(2, 1); dir[2] = eigenvectors(2, 2); dir.normalize(); float scale = 2/eigenvalues[2]; dir[0] = m_F*dir[0] + (1-m_F)*( (1-m_G)*dirOld[0] + scale*m_G*(tensor[0]*dirOld[0] + tensor[1]*dirOld[1] + tensor[2]*dirOld[2])); dir[1] = m_F*dir[1] + (1-m_F)*( (1-m_G)*dirOld[1] + scale*m_G*(tensor[1]*dirOld[0] + tensor[3]*dirOld[1] + tensor[4]*dirOld[2])); dir[2] = m_F*dir[2] + (1-m_F)*( (1-m_G)*dirOld[2] + scale*m_G*(tensor[2]*dirOld[0] + tensor[4]*dirOld[1] + tensor[5]*dirOld[2])); dir.normalize(); float angle = dot_product(dirOld, dir); if (angle<0) { dir *= -1; angle *= -1; } float r = m_StepSize/(2*std::asin(std::acos(angle)/2)); if (r void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, int threadId) { FiberPolyDataType poly = m_PolyDataContainer->GetElement(threadId); vtkSmartPointer points = vtkPoints::New(); vtkSmartPointer Cells = vtkCellArray::New(); typedef itk::DiffusionTensor3D TensorType; typedef ImageRegionConstIterator< InputImageType > InputIteratorType; typedef ImageRegionConstIterator< ItkUcharImgType > MaskIteratorType; typedef ImageRegionConstIterator< ItkFloatImgType > FloatIteratorType; typedef typename InputImageType::PixelType InputTensorType; InputIteratorType it(m_InputImage, outputRegionForThread ); MaskIteratorType mit(m_SeedImage, outputRegionForThread ); FloatIteratorType fit(m_FaImage, outputRegionForThread ); MaskIteratorType mit2(m_MaskImage, outputRegionForThread ); it.GoToBegin(); mit.GoToBegin(); mit2.GoToBegin(); fit.GoToBegin(); itk::Point worldPos; while( !it.IsAtEnd() ) { if (mit.Value()==0 || fit.Value() line = vtkSmartPointer::New(); std::vector< vtkIdType > pointIDs; typename InputImageType::IndexType index = it.GetIndex(); itk::ContinuousIndex start; unsigned int counter = 0; if (m_SeedsPerVoxel>1) { start[0] = index[0]+(double)(rand()%99-49)/100; start[1] = index[1]+(double)(rand()%99-49)/100; start[2] = index[2]+(double)(rand()%99-49)/100; } else { start[0] = index[0]; start[1] = index[1]; start[2] = index[2]; } // forward tracking float tractLength = FollowStreamline(start, 1, points, pointIDs); // add ids to line counter += pointIDs.size(); while (!pointIDs.empty()) { line->GetPointIds()->InsertNextId(pointIDs.back()); pointIDs.pop_back(); } // insert start point m_InputImage->TransformContinuousIndexToPhysicalPoint( start, worldPos ); line->GetPointIds()->InsertNextId(points->InsertNextPoint(worldPos.GetDataPointer())); // backward tracking tractLength += FollowStreamline(start, -1, points, pointIDs); counter += pointIDs.size(); if (tractLengthGetPointIds()->InsertNextId(pointIDs.at(i)); Cells->InsertNextCell(line); } ++mit; ++mit2; ++it; ++fit; } poly->SetPoints(points); poly->SetLines(Cells); std::cout << "Thread " << threadId << " finished tracking" << std::endl; } template< class TTensorPixelType, class TPDPixelType> vtkSmartPointer< vtkPolyData > StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2) { vtkSmartPointer vNewPolyData = vtkSmartPointer::New(); vtkSmartPointer vNewLines = poly1->GetLines(); vtkSmartPointer vNewPoints = poly1->GetPoints(); vtkSmartPointer vLines = poly2->GetLines(); vLines->InitTraversal(); for( int i=0; iGetNumberOfCells(); i++ ) { vtkIdType numPoints(0); vtkIdType* points(NULL); vLines->GetNextCell ( numPoints, points ); vtkSmartPointer container = vtkSmartPointer::New(); for( int j=0; jInsertNextPoint(poly2->GetPoint(points[j])); container->GetPointIds()->InsertNextId(id); } vNewLines->InsertNextCell(container); } // initialize polydata vNewPolyData->SetPoints(vNewPoints); vNewPolyData->SetLines(vNewLines); return vNewPolyData; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::AfterThreadedGenerateData() { MITK_INFO << "Generating polydata "; m_FiberPolyData = m_PolyDataContainer->GetElement(0); for (int i=1; iGetNumberOfThreads(); i++) { m_FiberPolyData = AddPolyData(m_FiberPolyData, m_PolyDataContainer->GetElement(i)); } MITK_INFO << "done"; } template< class TTensorPixelType, class TPDPixelType> void StreamlineTrackingFilter< TTensorPixelType, TPDPixelType> ::PrintSelf(std::ostream& os, Indent indent) const { } } #endif // __itkDiffusionQballPrincipleDirectionsImageFilter_txx diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index e1ba28e8cb..45a0e805e1 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,241 +1,242 @@ set(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures -> Tbss IODataStructures/TbssImages/mitkTbssImageSource.cpp IODataStructures/TbssImages/mitkTbssRoiImageSource.cpp IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageIOFactory.cpp IODataStructures/TbssImages/mitkTbssImage.cpp IODataStructures/TbssImages/mitkTbssRoiImage.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriterFactory.cpp IODataStructures/TbssImages/mitkTbssImporter.cpp # DataStructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp Rendering/mitkTbssImageMapper.cpp Rendering/mitkPlanarCircleMapper3D.cpp Rendering/mitkPlanarPolygonMapper3D.cpp Rendering/mitkConnectomicsNetworkMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/mitkTractAnalyzer.cpp # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.cpp Algorithms/Connectomics/mitkConnectomicsHistogramBase.cpp Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.cpp Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.cpp Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.cpp Algorithms/Connectomics/mitkConnectomicsHistogramCache.cpp Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.cpp # Tractography Tractography/GibbsTracking/mitkParticleGrid.cpp Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp Tractography/GibbsTracking/mitkEnergyComputer.cpp Tractography/GibbsTracking/mitkGibbsEnergyComputer.cpp Tractography/GibbsTracking/mitkFiberBuilder.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkTbssImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h Rendering/mitkPlanarCircleMapper3D.h Rendering/mitkPlanarPolygonMapper3D.h Rendering/mitkConnectomicsNetworkMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h Reconstruction/itkRegularizedIVIMReconstructionFilter.h Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/TbssImages/mitkTbssImporter.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h # Datastructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.h # Tractography Tractography/itkGibbsTrackingFilter.h Tractography/itkStochasticTractographyFilter.h Tractography/itkStreamlineTrackingFilter.h Tractography/GibbsTracking/mitkParticle.h Tractography/GibbsTracking/mitkParticleGrid.h Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h Tractography/GibbsTracking/mitkSimpSamp.h Tractography/GibbsTracking/mitkEnergyComputer.h Tractography/GibbsTracking/mitkGibbsEnergyComputer.h Tractography/GibbsTracking/mitkSphereInterpolator.h Tractography/GibbsTracking/mitkFiberBuilder.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkSkeletonizationFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h + Algorithms/itkMergeDiffusionImagesFilter.h # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.h Algorithms/Connectomics/mitkConnectomicsHistogramBase.h Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.h Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.h Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.h Algorithms/Connectomics/mitkConnectomicsHistogramCache.h Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.h ) set( TOOL_FILES ) if(WIN32) endif(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp ) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp index fbb9f9eb44..7ffc32adae 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp @@ -1,747 +1,747 @@ /*=================================================================== 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 // Qmitk #include "QmitkGibbsTrackingView.h" #include // Qt #include #include #include // MITK #include #include #include #include #include // ITK #include #include #include // MISC #include QmitkTrackingWorker::QmitkTrackingWorker(QmitkGibbsTrackingView* view) : m_View(view) { } void QmitkTrackingWorker::run() { m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); m_View->m_GlobalTracker->SetQBallImage(m_View->m_ItkQBallImage); m_View->m_GlobalTracker->SetTensorImage(m_View->m_ItkTensorImage); m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); m_View->m_GlobalTracker->SetStartTemperature((float)m_View->m_Controls->m_StartTempSlider->value()/100); m_View->m_GlobalTracker->SetEndTemperature((float)m_View->m_Controls->m_EndTempSlider->value()/10000); m_View->m_GlobalTracker->SetIterations(m_View->m_Iterations); m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); m_View->m_GlobalTracker->SetMinFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); m_View->m_GlobalTracker->SetCurvatureThreshold(cos((float)m_View->m_Controls->m_CurvatureThresholdSlider->value()*M_PI/180)); m_View->m_GlobalTracker->SetRandomSeed(m_View->m_Controls->m_RandomSeedSlider->value()); m_View->m_GlobalTracker->Update(); m_View->m_TrackingThread.quit(); } const std::string QmitkGibbsTrackingView::VIEW_ID = "org.mitk.views.gibbstracking"; QmitkGibbsTrackingView::QmitkGibbsTrackingView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) , m_ThreadIsRunning(false) , m_GlobalTracker(NULL) , m_QBallImage(NULL) , m_MaskImage(NULL) , m_ImageNode(NULL) , m_ItkQBallImage(NULL) , m_ItkTensorImage(NULL) , m_FiberBundleNode(NULL) , m_MaskImageNode(NULL) , m_TrackingWorker(this) , m_Iterations(10000000) , m_LastStep(0) { m_TrackingWorker.moveToThread(&m_TrackingThread); connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); connect(&m_TrackingThread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_TrackingTimer = new QTimer(this); } QmitkGibbsTrackingView::~QmitkGibbsTrackingView() { delete m_TrackingTimer; } // update tracking status and generate fiber bundle void QmitkGibbsTrackingView::TimerUpdate() { int currentStep = m_GlobalTracker->GetCurrentStep(); mitk::ProgressBar::GetInstance()->Progress(currentStep-m_LastStep); UpdateTrackingStatus(); GenerateFiberBundle(); m_LastStep = currentStep; } // tell global tractography filter to stop after current step void QmitkGibbsTrackingView::StopGibbsTracking() { if (m_GlobalTracker.IsNull()) return; //mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetSteps()-m_LastStep+1); m_GlobalTracker->SetAbortTracking(true); m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); } // update gui elements and generate fiber bundle after tracking is finished void QmitkGibbsTrackingView::AfterThread() { m_ThreadIsRunning = false; m_TrackingTimer->stop(); mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetSteps()-m_LastStep+1); UpdateGUI(); UpdateTrackingStatus(); if(m_Controls->m_ParticleWeightSlider->value()==0) { m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); } if(m_Controls->m_ParticleWidthSlider->value()==0) { m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); } if(m_Controls->m_ParticleLengthSlider->value()==0) { m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); } GenerateFiberBundle(); m_FiberBundleNode = NULL; } // start tracking timer and update gui elements before tracking is started void QmitkGibbsTrackingView::BeforeThread() { m_ThreadIsRunning = true; m_TrackingTime = QTime::currentTime(); m_ElapsedTime = 0; m_TrackingTimer->start(1000); m_LastStep = 0; UpdateGUI(); } // setup gui elements and signal/slot connections void QmitkGibbsTrackingView::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::QmitkGibbsTrackingViewControls; m_Controls->setupUi( parent ); AdvancedSettings(); connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); connect( m_Controls->m_IterationsSlider, SIGNAL(valueChanged(int)), this, SLOT(SetIterations(int)) ); connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); connect( m_Controls->m_CurvatureThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(SetCurvatureThreshold(int)) ); connect( m_Controls->m_RandomSeedSlider, SIGNAL(valueChanged(int)), this, SLOT(SetRandomSeed(int)) ); connect( m_Controls->m_OutputFileButton, SIGNAL(clicked()), this, SLOT(SetOutputFile()) ); } } void QmitkGibbsTrackingView::SetInExBalance(int value) { m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); } void QmitkGibbsTrackingView::SetFiberLength(int value) { m_Controls->m_FiberLengthLabel->setText(QString::number(value)+"mm"); } void QmitkGibbsTrackingView::SetRandomSeed(int value) { if (value>=0) m_Controls->m_RandomSeedLabel->setText(QString::number(value)); else m_Controls->m_RandomSeedLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleWeight(int value) { if (value>0) m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); else m_Controls->m_ParticleWeightLabel->setText("auto"); } void QmitkGibbsTrackingView::SetStartTemp(int value) { m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); } void QmitkGibbsTrackingView::SetEndTemp(int value) { m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); } void QmitkGibbsTrackingView::SetParticleWidth(int value) { if (value>0) m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleWidthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleLength(int value) { if (value>0) m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); else m_Controls->m_ParticleLengthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetCurvatureThreshold(int value) { m_Controls->m_CurvatureThresholdLabel->setText(QString::number(value)+"°"); } void QmitkGibbsTrackingView::SetIterations(int value) { switch(value) { case 0: m_Controls->m_IterationsLabel->setText("Iterations: 1x10^4"); m_Iterations = 10000; break; case 1: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); m_Iterations = 50000; break; case 2: m_Controls->m_IterationsLabel->setText("Iterations: 1x10^5"); m_Iterations = 100000; break; case 3: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); m_Iterations = 500000; break; case 4: m_Controls->m_IterationsLabel->setText("Iterations: 1x10^6"); m_Iterations = 1000000; break; case 5: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); m_Iterations = 5000000; break; case 6: m_Controls->m_IterationsLabel->setText("Iterations: 1x10^7"); m_Iterations = 10000000; break; case 7: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); m_Iterations = 50000000; break; case 8: m_Controls->m_IterationsLabel->setText("Iterations: 1x10^8"); m_Iterations = 100000000; break; case 9: m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); m_Iterations = 500000000; break; case 10: m_Controls->m_IterationsLabel->setText("Iterations: 1x10^9"); m_Iterations = 1000000000; break; } } void QmitkGibbsTrackingView::StdMultiWidgetAvailable(QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkGibbsTrackingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } // called if datamanager selection changes void QmitkGibbsTrackingView::OnSelectionChanged( std::vector nodes ) { if (m_ThreadIsRunning) return; m_ImageNode = NULL; m_MaskImageNode = NULL; // iterate all selected objects for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast(node->GetData()) ) m_ImageNode = node; else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) m_ImageNode = node; else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) m_MaskImageNode = node; } } UpdateGUI(); } // update gui elements displaying trackings status void QmitkGibbsTrackingView::UpdateTrackingStatus() { if (m_GlobalTracker.IsNull()) return; m_ElapsedTime += m_TrackingTime.elapsed()/1000; m_TrackingTime.restart(); unsigned long hours = m_ElapsedTime/3600; unsigned long minutes = (m_ElapsedTime%3600)/60; unsigned long seconds = m_ElapsedTime%60; m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); m_Controls->m_CurrentStepLabel->setText( QString::number(100*(float)(m_GlobalTracker->GetCurrentStep()-1)/m_GlobalTracker->GetSteps())+"%" ); m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); } // update gui elements (enable/disable elements and set tooltips) void QmitkGibbsTrackingView::UpdateGUI() { if (m_ImageNode.IsNotNull()) m_Controls->m_QballImageLabel->setText(m_ImageNode->GetName().c_str()); else m_Controls->m_QballImageLabel->setText("-"); if (m_MaskImageNode.IsNotNull()) m_Controls->m_MaskImageLabel->setText(m_MaskImageNode->GetName().c_str()); else m_Controls->m_MaskImageLabel->setText("-"); if (!m_ThreadIsRunning && m_ImageNode.IsNotNull()) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(true); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_IterationsSlider->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); m_Controls->m_TrackingStop->setToolTip(""); } else if (!m_ThreadIsRunning) { m_Controls->m_TrackingStop->setEnabled(false); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(true); m_Controls->m_IterationsSlider->setEnabled(true); m_Controls->m_AdvancedFrame->setEnabled(true); m_Controls->m_TrackingStop->setText("Stop Tractography"); m_Controls->m_TrackingStart->setToolTip("No Q-Ball image selected."); m_Controls->m_TrackingStop->setToolTip(""); } else { m_Controls->m_TrackingStop->setEnabled(true); m_Controls->m_TrackingStart->setEnabled(false); m_Controls->m_LoadTrackingParameters->setEnabled(false); m_Controls->m_IterationsSlider->setEnabled(false); m_Controls->m_AdvancedFrame->setEnabled(false); m_Controls->m_AdvancedFrame->setVisible(false); m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); } } // show/hide advanced settings frame void QmitkGibbsTrackingView::AdvancedSettings() { m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); } // set mask image data node void QmitkGibbsTrackingView::SetMask() { std::vector nodes = GetDataManagerSelection(); if (nodes.empty()) { m_MaskImageNode = NULL; m_Controls->m_MaskImageLabel->setText("-"); return; } for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if (node.IsNotNull() && dynamic_cast(node->GetData())) { m_MaskImageNode = node; m_Controls->m_MaskImageLabel->setText(node->GetName().c_str()); return; } } } // cast image to float template void QmitkGibbsTrackingView::CastToFloat(InputImageType* image, mitk::Image::Pointer outImage) { typedef itk::CastImageFilter ItkCastFilter; typename ItkCastFilter::Pointer itkCaster = ItkCastFilter::New(); itkCaster->SetInput(image); itkCaster->Update(); outImage->InitializeByItk(itkCaster->GetOutput()); outImage->SetVolume(itkCaster->GetOutput()->GetBufferPointer()); } // check for mask and qbi and start tracking thread void QmitkGibbsTrackingView::StartGibbsTracking() { if(m_ThreadIsRunning) { MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; return; } if (m_ImageNode.IsNull()) { QMessageBox::information( NULL, "Warning", "Please load and select a qball image before starting image processing."); return; } if (dynamic_cast(m_ImageNode->GetData())) m_QBallImage = dynamic_cast(m_ImageNode->GetData()); else if (dynamic_cast(m_ImageNode->GetData())) m_TensorImage = dynamic_cast(m_ImageNode->GetData()); if (m_QBallImage.IsNull() && m_TensorImage.IsNull()) return; // cast qbi to itk m_ItkTensorImage = NULL; m_ItkQBallImage = NULL; m_MaskImage = NULL; if (m_QBallImage.IsNotNull()) { m_ItkQBallImage = ItkQBallImgType::New(); mitk::CastToItkImage(m_QBallImage, m_ItkQBallImage); } else { m_ItkTensorImage = ItkTensorImage::New(); mitk::CastToItkImage(m_TensorImage, m_ItkTensorImage); } // mask image found? // catch exceptions thrown by the itkAccess macros try{ if(m_MaskImageNode.IsNotNull()) { if (dynamic_cast(m_MaskImageNode->GetData())) mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); } } catch(...){}; unsigned int steps = m_Iterations/10000; if (steps<10) steps = 10; m_LastStep = 1; mitk::ProgressBar::GetInstance()->AddStepsToDo(steps); // start worker thread m_TrackingThread.start(QThread::LowestPriority); } // generate mitkFiberBundle from tracking filter output void QmitkGibbsTrackingView::GenerateFiberBundle() { if (m_GlobalTracker.IsNull() || (!(m_Controls->m_VisualizationCheckbox->isChecked() || m_Controls->m_VisualizeOnceButton->isChecked()) && m_ThreadIsRunning)) return; if (m_Controls->m_VisualizeOnceButton->isChecked()) m_Controls->m_VisualizeOnceButton->setChecked(false); vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); if ( m_GlobalTracker->GetNumAcceptedFibers()==0 ) return; m_FiberBundle = mitk::FiberBundleX::New(fiberBundle); if (m_FiberBundleNode.IsNotNull()){ GetDefaultDataStorage()->Remove(m_FiberBundleNode); m_FiberBundleNode = 0; } m_FiberBundleNode = mitk::DataNode::New(); m_FiberBundleNode->SetData(m_FiberBundle); QString name(m_ImageNode->GetName().c_str()); - name += "_FiberBundle"; + name += "_Gibbs"; m_FiberBundleNode->SetName(name.toStdString()); m_FiberBundleNode->SetVisibility(true); if (!m_OutputFileName.isEmpty()) { QString filename = m_OutputFileName; mitk::FiberBundleXWriter::Pointer writer = mitk::FiberBundleXWriter::New(); writer->SetFileName(filename.toStdString()); writer->SetInputFiberBundleX(m_FiberBundle.GetPointer()); try { writer->Update(); QMessageBox::information(NULL, "Fiber bundle saved to", filename); } catch (itk::ExceptionObject &ex) { QMessageBox::information(NULL, "Fiber bundle could not be saved", QString("%1\n%2\n%3\n%4\n%5\n%6").arg(ex.GetNameOfClass()).arg(ex.GetFile()).arg(ex.GetLine()).arg(ex.GetLocation()).arg(ex.what()).arg(ex.GetDescription())); if(m_ImageNode.IsNull()) GetDataStorage()->Add(m_FiberBundleNode); else GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); } } else { if(m_ImageNode.IsNull()) GetDataStorage()->Add(m_FiberBundleNode); else GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); } } void QmitkGibbsTrackingView::SetOutputFile() { // SELECT FOLDER DIALOG m_OutputFileName = QFileDialog::getSaveFileName(0, tr("Set file name"), QDir::currentPath()+"/FiberBundle.fib", tr("Fiber Bundle (*.fib)") ); if (m_OutputFileName.isEmpty()) m_Controls->m_OutputFileLabel->setText("N/A"); else m_Controls->m_OutputFileLabel->setText(m_OutputFileName); } // save current tracking paramters as xml file (.gtp) void QmitkGibbsTrackingView::SaveTrackingParameters() { TiXmlDocument documentXML; TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); documentXML.LinkEndChild( declXML ); TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); mainXML->SetAttribute("file_version", "0.1"); documentXML.LinkEndChild(mainXML); TiXmlElement* paramXML = new TiXmlElement("parameter_set"); paramXML->SetAttribute("iterations", QString::number(m_Iterations).toStdString()); paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); paramXML->SetAttribute("curvature_threshold", QString::number(m_Controls->m_CurvatureThresholdSlider->value()).toStdString()); mainXML->LinkEndChild(paramXML); QString filename = QFileDialog::getSaveFileName( 0, tr("Save Parameters"), QDir::currentPath()+"/param.gtp", tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; if(!filename.endsWith(".gtp")) filename += ".gtp"; documentXML.SaveFile( filename.toStdString() ); } void QmitkGibbsTrackingView::UpdateIteraionsGUI(unsigned long iterations) { switch(iterations) { case 10000: m_Controls->m_IterationsSlider->setValue(0); m_Controls->m_IterationsLabel->setText("Iterations: 10^4"); break; case 50000: m_Controls->m_IterationsSlider->setValue(1); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); break; case 100000: m_Controls->m_IterationsSlider->setValue(2); m_Controls->m_IterationsLabel->setText("Iterations: 10^5"); break; case 500000: m_Controls->m_IterationsSlider->setValue(3); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); break; case 1000000: m_Controls->m_IterationsSlider->setValue(4); m_Controls->m_IterationsLabel->setText("Iterations: 10^6"); break; case 5000000: m_Controls->m_IterationsSlider->setValue(5); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); break; case 10000000: m_Controls->m_IterationsSlider->setValue(6); m_Controls->m_IterationsLabel->setText("Iterations: 10^7"); break; case 50000000: m_Controls->m_IterationsSlider->setValue(7); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); break; case 100000000: m_Controls->m_IterationsSlider->setValue(8); m_Controls->m_IterationsLabel->setText("Iterations: 10^8"); break; case 500000000: m_Controls->m_IterationsSlider->setValue(9); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); break; case 1000000000: m_Controls->m_IterationsSlider->setValue(10); m_Controls->m_IterationsLabel->setText("Iterations: 10^9"); break; case 5000000000: m_Controls->m_IterationsSlider->setValue(11); m_Controls->m_IterationsLabel->setText("Iterations: 5x10^9"); break; } } // load current tracking paramters from xml file (.gtp) void QmitkGibbsTrackingView::LoadTrackingParameters() { QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); if(filename.isEmpty() || filename.isNull()) return; TiXmlDocument doc( filename.toStdString() ); doc.LoadFile(); TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("parameter_set").Element(); QString iterations(pElem->Attribute("iterations")); m_Iterations = iterations.toULong(); UpdateIteraionsGUI(m_Iterations); QString particleLength(pElem->Attribute("particle_length")); float pLength = particleLength.toFloat(); QString particleWidth(pElem->Attribute("particle_width")); float pWidth = particleWidth.toFloat(); if (pLength==0) m_Controls->m_ParticleLengthLabel->setText("auto"); else m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); if (pWidth==0) m_Controls->m_ParticleWidthLabel->setText("auto"); else m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); m_Controls->m_ParticleLengthSlider->setValue(pLength*10); QString partWeight(pElem->Attribute("particle_weight")); m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); m_Controls->m_ParticleWeightLabel->setText(partWeight); QString startTemp(pElem->Attribute("temp_start")); m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); m_Controls->m_StartTempLabel->setText(startTemp); QString endTemp(pElem->Attribute("temp_end")); m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); m_Controls->m_EndTempLabel->setText(endTemp); QString inExBalance(pElem->Attribute("inexbalance")); m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); m_Controls->m_InExBalanceLabel->setText(inExBalance); QString fiberLength(pElem->Attribute("fiber_length")); m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); m_Controls->m_FiberLengthLabel->setText(fiberLength+"mm"); QString curvThres(pElem->Attribute("curvature_threshold")); m_Controls->m_CurvatureThresholdSlider->setValue(curvThres.toInt()); m_Controls->m_CurvatureThresholdLabel->setText(curvThres+"°"); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp index 5b8a98a47d..34a8735298 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.cpp @@ -1,580 +1,643 @@ /*=================================================================== 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. ===================================================================*/ //#define MBILOG_ENABLE_DEBUG #include "QmitkPreprocessingView.h" #include "mitkDiffusionImagingConfigure.h" // qt includes #include // itk includes #include "itkTimeProbe.h" #include "itkB0ImageExtractionImageFilter.h" #include "itkB0ImageExtractionToSeparateImageFilter.h" #include "itkBrainMaskExtractionImageFilter.h" #include "itkCastImageFilter.h" #include "itkVectorContainer.h" #include +#include // mitk includes #include "QmitkDataStorageComboBox.h" #include "QmitkStdMultiWidget.h" #include "mitkProgressBar.h" #include "mitkStatusBar.h" #include "mitkNodePredicateDataType.h" #include "mitkProperties.h" #include "mitkVtkResliceInterpolationProperty.h" #include "mitkLookupTable.h" #include "mitkLookupTableProperty.h" #include "mitkTransferFunction.h" #include "mitkTransferFunctionProperty.h" #include "mitkDataNodeObject.h" #include "mitkOdfNormalizationMethodProperty.h" #include "mitkOdfScaleByProperty.h" #include #include #include const std::string QmitkPreprocessingView::VIEW_ID = "org.mitk.views.preprocessing"; #define DI_INFO MITK_INFO("DiffusionImaging") typedef float TTensorPixelType; QmitkPreprocessingView::QmitkPreprocessingView() : QmitkFunctionality(), m_Controls(NULL), m_MultiWidget(NULL), m_DiffusionImage(NULL) { } QmitkPreprocessingView::QmitkPreprocessingView(const QmitkPreprocessingView& other) { Q_UNUSED(other) throw std::runtime_error("Copy constructor not implemented"); } QmitkPreprocessingView::~QmitkPreprocessingView() { } void QmitkPreprocessingView::CreateQtPartControl(QWidget *parent) { if (!m_Controls) { // create GUI widgets m_Controls = new Ui::QmitkPreprocessingViewControls; m_Controls->setupUi(parent); this->CreateConnections(); m_Controls->m_MeasurementFrameTable->horizontalHeader()->setResizeMode(QHeaderView::Stretch); m_Controls->m_MeasurementFrameTable->verticalHeader()->setResizeMode(QHeaderView::Stretch); } } void QmitkPreprocessingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkPreprocessingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkPreprocessingView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_ButtonAverageGradients), SIGNAL(clicked()), this, SLOT(AverageGradients()) ); connect( (QObject*)(m_Controls->m_ButtonExtractB0), SIGNAL(clicked()), this, SLOT(ExtractB0()) ); connect( (QObject*)(m_Controls->m_ButtonBrainMask), SIGNAL(clicked()), this, SLOT(BrainMask()) ); connect( (QObject*)(m_Controls->m_ModifyMeasurementFrame), SIGNAL(clicked()), this, SLOT(DoApplyMesurementFrame()) ); connect( (QObject*)(m_Controls->m_ReduceGradientsButton), SIGNAL(clicked()), this, SLOT(DoReduceGradientDirections()) ); connect( (QObject*)(m_Controls->m_ShowGradientsButton), SIGNAL(clicked()), this, SLOT(DoShowGradientDirections()) ); connect( (QObject*)(m_Controls->m_MirrorGradientToHalfSphereButton), SIGNAL(clicked()), this, SLOT(DoHalfSphereGradientDirections()) ); + connect( (QObject*)(m_Controls->m_MergeDwisButton), SIGNAL(clicked()), this, SLOT(MergeDwis()) ); + } } void QmitkPreprocessingView::OnSelectionChanged( std::vector nodes ) { bool foundDwiVolume = false; m_DiffusionImage = NULL; m_SelectedDiffusionNodes.clear(); // iterate selection for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast*>(node->GetData()) ) { foundDwiVolume = true; m_DiffusionImage = dynamic_cast*>(node->GetData()); m_Controls->m_DiffusionImageLabel->setText(node->GetName().c_str()); m_SelectedDiffusionNodes.push_back(node); } } m_Controls->m_ButtonBrainMask->setEnabled(foundDwiVolume); m_Controls->m_ButtonAverageGradients->setEnabled(foundDwiVolume); m_Controls->m_ButtonExtractB0->setEnabled(foundDwiVolume); m_Controls->m_CheckExtractAll->setEnabled(foundDwiVolume); m_Controls->m_ModifyMeasurementFrame->setEnabled(foundDwiVolume); m_Controls->m_MeasurementFrameTable->setEnabled(foundDwiVolume); m_Controls->m_ReduceGradientsButton->setEnabled(foundDwiVolume); m_Controls->m_ShowGradientsButton->setEnabled(foundDwiVolume); m_Controls->m_MirrorGradientToHalfSphereButton->setEnabled(foundDwiVolume); + m_Controls->m_MergeDwisButton->setEnabled(foundDwiVolume); foreach(QCheckBox * box, m_ReduceGradientCheckboxes) { m_Controls->m_ReductionFrame->layout()->removeWidget(box); delete box; } foreach(QSpinBox * box, m_ReduceGradientSpinboxes) { m_Controls->m_ReductionFrame->layout()->removeWidget(box); delete box; } m_ReduceGradientCheckboxes.clear(); m_ReduceGradientSpinboxes.clear(); if (foundDwiVolume) { vnl_matrix_fixed< double, 3, 3 > mf = m_DiffusionImage->GetMeasurementFrame(); for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); delete item; item = new QTableWidgetItem(); item->setTextAlignment(Qt::AlignCenter | Qt::AlignVCenter); item->setText(QString::number(mf.get(r,c))); m_Controls->m_MeasurementFrameTable->setItem(r,c,item); } typedef mitk::DiffusionImage::BValueMap BValueMap; typedef mitk::DiffusionImage::BValueMap::iterator BValueMapIterator; BValueMap bValMap = m_DiffusionImage->GetB_ValueMap(); BValueMapIterator it = bValMap.begin(); m_Controls->m_BvalueTable->clear(); m_Controls->m_BvalueTable->setRowCount(bValMap.size() ); QStringList headerList; headerList << "b-Value" << "Number of gradients"; m_Controls->m_BvalueTable->setHorizontalHeaderLabels(headerList); QCheckBox* checkBox; QSpinBox* spinBox; int i = 0 ; for(;it != bValMap.end(); it++) { m_Controls->m_BvalueTable->setItem(i,0,new QTableWidgetItem(QString::number(it->first))); m_Controls->m_BvalueTable->setItem(i,1,new QTableWidgetItem(QString::number(it->second.size()))); // Reduce Gradients GUI adaption if(it->first != 0 && bValMap.size() > 1){ checkBox = new QCheckBox(QString::number(it->first) + " with " + QString::number(it->second.size()) + " directions"); checkBox->setEnabled(true); checkBox->setChecked(true); checkBox->setCheckable(true); m_ReduceGradientCheckboxes.push_back(checkBox); m_Controls->m_ReductionFrame->layout()->addWidget(checkBox); spinBox = new QSpinBox(); spinBox->setValue(std::ceil((float)it->second.size()/2)); spinBox->setMaximum(it->second.size()); spinBox->setMinimum(0); m_ReduceGradientSpinboxes.push_back(spinBox); m_Controls->m_ReductionFrame->layout()->addWidget(spinBox); } i++; } } else { for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); delete item; item = new QTableWidgetItem(); m_Controls->m_MeasurementFrameTable->setItem(r,c,item); } m_Controls->m_BvalueTable->clear(); m_Controls->m_BvalueTable->setRowCount(1); QStringList headerList; headerList << "b-Value" << "Number of gradients"; m_Controls->m_BvalueTable->setHorizontalHeaderLabels(headerList); m_Controls->m_BvalueTable->setItem(0,0,new QTableWidgetItem("-")); m_Controls->m_BvalueTable->setItem(0,1,new QTableWidgetItem("-")); m_Controls->m_DiffusionImageLabel->setText("-"); } } void QmitkPreprocessingView::Activated() { QmitkFunctionality::Activated(); } void QmitkPreprocessingView::Deactivated() { QmitkFunctionality::Deactivated(); } void QmitkPreprocessingView::DoHalfSphereGradientDirections() { GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetDirections(); for (int j=0; jSize(); j++) if (gradientContainer->at(j)[0]<0) gradientContainer->at(j) = -gradientContainer->at(j); m_DiffusionImage->SetDirections(gradientContainer); } void QmitkPreprocessingView::DoApplyMesurementFrame() { if (m_DiffusionImage.IsNull()) return; vnl_matrix_fixed< double, 3, 3 > mf; for (int r=0; r<3; r++) for (int c=0; c<3; c++) { QTableWidgetItem* item = m_Controls->m_MeasurementFrameTable->item(r,c); if (!item) return; mf[r][c] = item->text().toDouble(); } m_DiffusionImage->SetMeasurementFrame(mf); } void QmitkPreprocessingView::DoShowGradientDirections() { if (m_DiffusionImage.IsNull()) return; int maxIndex = 0; int maxSize = m_DiffusionImage->GetDimension(0); if (maxSizeGetDimension(1)) { maxSize = m_DiffusionImage->GetDimension(1); maxIndex = 1; } if (maxSizeGetDimension(2)) { maxSize = m_DiffusionImage->GetDimension(2); maxIndex = 2; } mitk::Point3D origin = m_DiffusionImage->GetGeometry()->GetOrigin(); mitk::PointSet::Pointer originSet = mitk::PointSet::New(); typedef mitk::DiffusionImage::BValueMap BValueMap; typedef mitk::DiffusionImage::BValueMap::iterator BValueMapIterator; BValueMap bValMap = m_DiffusionImage->GetB_ValueMap(); GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetDirectionsWithMeasurementFrame(); mitk::Geometry3D::Pointer geometry = m_DiffusionImage->GetGeometry(); int shellCount = 1; for(BValueMapIterator it = bValMap.begin(); it!=bValMap.end(); ++it) { mitk::PointSet::Pointer pointset = mitk::PointSet::New(); for (int j=0; jsecond.size(); j++) { mitk::Point3D ip; vnl_vector_fixed< double, 3 > v = gradientContainer->at(it->second[j]); if (v.magnitude()>mitk::eps) { ip[0] = v[0]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[0]-0.5*geometry->GetSpacing()[0] + geometry->GetSpacing()[0]*m_DiffusionImage->GetDimension(0)/2; ip[1] = v[1]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[1]-0.5*geometry->GetSpacing()[1] + geometry->GetSpacing()[1]*m_DiffusionImage->GetDimension(1)/2; ip[2] = v[2]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[2]-0.5*geometry->GetSpacing()[2] + geometry->GetSpacing()[2]*m_DiffusionImage->GetDimension(2)/2; pointset->InsertPoint(j, ip); } else if (originSet->IsEmpty()) { ip[0] = v[0]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[0]-0.5*geometry->GetSpacing()[0] + geometry->GetSpacing()[0]*m_DiffusionImage->GetDimension(0)/2; ip[1] = v[1]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[1]-0.5*geometry->GetSpacing()[1] + geometry->GetSpacing()[1]*m_DiffusionImage->GetDimension(1)/2; ip[2] = v[2]*maxSize*geometry->GetSpacing()[maxIndex]/2 + origin[2]-0.5*geometry->GetSpacing()[2] + geometry->GetSpacing()[2]*m_DiffusionImage->GetDimension(2)/2; originSet->InsertPoint(j, ip); } } if (it->firstSetData(pointset); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); name += "_Shell_"; name += QString::number(it->first); node->SetName(name.toStdString().c_str()); node->SetProperty("pointsize", mitk::FloatProperty::New((float)maxSize/50)); int b0 = shellCount%2; int b1 = 0; int b2 = 0; if (shellCount>4) b2 = 1; if (shellCount%4 >= 2) b1 = 1; node->SetProperty("color", mitk::ColorProperty::New(b2, b1, b0)); GetDefaultDataStorage()->Add(node, m_SelectedDiffusionNodes.front()); shellCount++; } // add origin to datastorage mitk::DataNode::Pointer node = mitk::DataNode::New(); node->SetData(originSet); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); name += "_Origin"; node->SetName(name.toStdString().c_str()); node->SetProperty("pointsize", mitk::FloatProperty::New((float)maxSize/50)); node->SetProperty("color", mitk::ColorProperty::New(1,1,1)); GetDefaultDataStorage()->Add(node, m_SelectedDiffusionNodes.front()); } void QmitkPreprocessingView::DoReduceGradientDirections() { if (m_DiffusionImage.IsNull()) return; typedef mitk::DiffusionImage DiffusionImageType; typedef itk::ElectrostaticRepulsionDiffusionGradientReductionFilter FilterType; typedef DiffusionImageType::BValueMap BValueMap; // GetShellSelection from GUI BValueMap shellSlectionMap; BValueMap originalShellMap = m_DiffusionImage->GetB_ValueMap(); std::vector newNumGradientDirections; int shellCounter = 0; foreach(QCheckBox * box , m_ReduceGradientCheckboxes) { if(box->isChecked()) { double BValue = (box->text().split(' ')).at(0).toDouble(); shellSlectionMap[BValue] = originalShellMap[BValue]; newNumGradientDirections.push_back(m_ReduceGradientSpinboxes.at(shellCounter)->value()); } shellCounter++; } if (newNumGradientDirections.empty()) return; GradientDirectionContainerType::Pointer gradientContainer = m_DiffusionImage->GetDirections(); FilterType::Pointer filter = FilterType::New(); filter->SetInput(m_DiffusionImage->GetVectorImage()); filter->SetOriginalGradientDirections(gradientContainer); filter->SetNumGradientDirections(newNumGradientDirections); filter->SetOriginalBValueMap(originalShellMap); filter->SetShellSelectionBValueMap(shellSlectionMap); filter->Update(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage( filter->GetOutput() ); image->SetB_Value(m_DiffusionImage->GetB_Value()); image->SetDirections(filter->GetGradientDirections()); image->SetMeasurementFrame(m_DiffusionImage->GetMeasurementFrame()); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); foreach(QSpinBox* box, m_ReduceGradientSpinboxes) { name += "_"; name += QString::number(box->value()); } imageNode->SetName(name.toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); } +void QmitkPreprocessingView::MergeDwis() +{ + typedef mitk::DiffusionImage DiffusionImageType; + typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; + + if (m_SelectedDiffusionNodes.size()<2) + return; + + typedef itk::VectorImage DwiImageType; + typedef typename DwiImageType::PixelType DwiPixelType; + typedef typename DwiImageType::RegionType DwiRegionType; + typedef typename std::vector< typename DwiImageType::Pointer > DwiImageContainerType; + + typedef typename std::vector< GradientContainerType::Pointer > GradientListContainerType; + + DwiImageContainerType imageContainer; + GradientListContainerType gradientListContainer; + std::vector< double > bValueContainer; + + for (int i=0; i* >( m_SelectedDiffusionNodes.at(i)->GetData() ); + if ( dwi.IsNotNull() ) + { + imageContainer.push_back(dwi->GetVectorImage()); + gradientListContainer.push_back(dwi->GetDirectionsWithMeasurementFrame()); + bValueContainer.push_back(dwi->GetB_Value()); + } + } + + typedef itk::MergeDiffusionImagesFilter FilterType; + FilterType::Pointer filter = FilterType::New(); + filter->SetImageVolumes(imageContainer); + filter->SetGradientLists(gradientListContainer); + filter->SetBValues(bValueContainer); + filter->Update(); + + vnl_matrix_fixed< double, 3, 3 > mf; mf.set_identity(); + DiffusionImageType::Pointer image = DiffusionImageType::New(); + image->SetVectorImage( filter->GetOutput() ); + image->SetB_Value(filter->GetBValue()); + image->SetDirections(filter->GetOutputGradients()); + image->SetMeasurementFrame(mf); + image->InitializeFromVectorImage(); + + mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); + imageNode->SetData( image ); + QString name = m_SelectedDiffusionNodes.front()->GetName().c_str(); + + for (int i=0; iSetName(name.toStdString().c_str()); + GetDefaultDataStorage()->Add(imageNode); +} + void QmitkPreprocessingView::ExtractB0() { typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; int nrFiles = m_SelectedDiffusionNodes.size(); if (!nrFiles) return; // call the extraction withou averaging if the check-box is checked if( this->m_Controls->m_CheckExtractAll->isChecked() ) { DoExtractBOWithoutAveraging(); return; } mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes.begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes.end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { DiffusionImageType* vols = static_cast( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); filter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( filter->GetOutput() ); mitkImage->SetVolume( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_B0")); GetDefaultDataStorage()->Add(node); ++itemiter; } } void QmitkPreprocessingView::DoExtractBOWithoutAveraging() { // typedefs typedef mitk::DiffusionImage DiffusionImageType; typedef DiffusionImageType::GradientDirectionContainerType GradientContainerType; typedef itk::B0ImageExtractionToSeparateImageFilter< short, short> FilterType; // check number of selected objects, return if empty int nrFiles = m_SelectedDiffusionNodes.size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes.begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes.end() ); while ( itemiter != itemiterend ) // for all items { DiffusionImageType* vols = static_cast( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); filter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( filter->GetOutput() ); mitkImage->SetImportChannel( filter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_B0_ALL")); GetDefaultDataStorage()->Add(node); ++itemiter; } } void QmitkPreprocessingView::AverageGradients() { int nrFiles = m_SelectedDiffusionNodes.size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes.begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes.end() ); std::vector nodes; while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); vols->AverageRedundantGradients(m_Controls->m_Blur->value()); ++itemiter; } } void QmitkPreprocessingView::BrainMask() { int nrFiles = m_SelectedDiffusionNodes.size(); if (!nrFiles) return; mitk::DataStorage::SetOfObjects::const_iterator itemiter( m_SelectedDiffusionNodes.begin() ); mitk::DataStorage::SetOfObjects::const_iterator itemiterend( m_SelectedDiffusionNodes.end() ); while ( itemiter != itemiterend ) // for all items { mitk::DiffusionImage* vols = static_cast*>( (*itemiter)->GetData()); std::string nodename; (*itemiter)->GetStringProperty("name", nodename); // Extract image using found index typedef itk::B0ImageExtractionImageFilter FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(vols->GetVectorImage()); filter->SetDirections(vols->GetDirections()); typedef itk::CastImageFilter, itk::Image > CastFilterType; CastFilterType::Pointer castfilter = CastFilterType::New(); castfilter->SetInput(filter->GetOutput()); typedef itk::BrainMaskExtractionImageFilter MaskFilterType; MaskFilterType::Pointer maskfilter = MaskFilterType::New(); maskfilter->SetInput(castfilter->GetOutput()); maskfilter->Update(); mitk::Image::Pointer mitkImage = mitk::Image::New(); mitkImage->InitializeByItk( maskfilter->GetOutput() ); mitkImage->SetVolume( maskfilter->GetOutput()->GetBufferPointer() ); mitk::DataNode::Pointer node=mitk::DataNode::New(); node->SetData( mitkImage ); node->SetProperty( "name", mitk::StringProperty::New(nodename + "_Mask")); GetDefaultDataStorage()->Add(node); ++itemiter; } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h index fa7eb63ff2..f6c3b6eea6 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingView.h @@ -1,110 +1,110 @@ /*=================================================================== 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 _QMITKPREPROCESSINGVIEW_H_INCLUDED #define _QMITKPREPROCESSINGVIEW_H_INCLUDED #include #include #include "ui_QmitkPreprocessingViewControls.h" #include "mitkDiffusionImage.h" typedef short DiffusionPixelType; struct PrpSelListener; /*! * \ingroup org_mitk_gui_qt_preprocessing_internal * * \brief QmitkPreprocessingView * * Document your class here. * * \sa QmitkFunctionality */ class QmitkPreprocessingView : public QmitkFunctionality { friend struct PrpSelListener; // this is needed for all Qt objects that should have a MOC object (everything that derives from QObject) Q_OBJECT public: static const std::string VIEW_ID; typedef vnl_vector_fixed< double, 3 > GradientDirectionType; typedef itk::VectorContainer< unsigned int, GradientDirectionType > GradientDirectionContainerType; QmitkPreprocessingView(); QmitkPreprocessingView(const QmitkPreprocessingView& other); virtual ~QmitkPreprocessingView(); virtual void CreateQtPartControl(QWidget *parent); /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); /// \brief Called when the functionality is activated virtual void Activated(); virtual void Deactivated(); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); static const int nrconvkernels; protected slots: void AverageGradients(); void ExtractB0(); void BrainMask(); - + void MergeDwis(); void DoApplyMesurementFrame(); void DoReduceGradientDirections(); void DoShowGradientDirections(); void DoHalfSphereGradientDirections(); protected: /** Called by ExtractB0 if check-box activated, extracts all b0 images without averaging */ void DoExtractBOWithoutAveraging(); /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); Ui::QmitkPreprocessingViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; void SetDefaultNodeProperties(mitk::DataNode::Pointer node, std::string name); mitk::DiffusionImage::Pointer m_DiffusionImage; std::vector< mitk::DataNode::Pointer > m_SelectedDiffusionNodes; QList m_ReduceGradientCheckboxes; QList m_ReduceGradientSpinboxes; }; #endif // _QMITKPREPROCESSINGVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui index bae365a5fd..3ae52c2eb5 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkPreprocessingViewControls.ui @@ -1,549 +1,577 @@ QmitkPreprocessingViewControls 0 0 389 992 0 0 false QmitkPreprocessingViewControls true Data Diffusion Image: - 0 0 Info 0 0 Qt::ScrollBarAsNeeded Qt::ScrollBarAlwaysOff true 100 true false true b-Value Number of gradients Qt::Horizontal 40 20 false Generate pointset displaying the gradient vectors (applied measurement frame). Show gradients Reduce size 0 70 Multiple acquistions of one gradient direction can be averaged. Due to rounding errors, similar gradients often differ in the last decimal positions. The Merge radius allows to average them anyway by taking into account all directions within a certain radius. true QFrame::NoFrame QFrame::Raised 0 Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. 6 2.000000000000000 0.000100000000000 0.001000000000000 Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Accumulates the information that was acquired with multiple repetitions for one gradient. Vectors do not have to be precisely equal in order to be merged, if a "Merge radius" > 0 is configured. Merge radius false Average redundant gradients false Mirror all gradients around one axis. Mirror gradients to half sphere Qt::Horizontal QFrame::NoFrame QFrame::Raised 9 9 Specify desired number of gradients per shell. false Retain only the specified number of gradient directions and according image volumes. The retained directions are spread equally over the half sphere. Reduce number of gradients Non diffusion weighted image 0 30 Average and extract all images that were acquired without diffusion weighting. true false Extract B0 Create a 3D+t data set containing all b0 images as timesteps Extract all B0 without averaging + + + + Merge selected images + + + + + + false + + + + + + + + + + + + Merge + + + + + + Brain mask false Estimate binary brain mask 0 0 Measurment frame Qt::Horizontal 40 20 false 0 0 10 10 IBeamCursor true Qt::ScrollBarAlwaysOff Qt::ScrollBarAlwaysOff true false false true true false true true New Row New Row New Row New Column New Column New Column false Apply new mesurement frame Qt::Vertical 20 40