diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp index 954a726f06..930a69efee 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp @@ -1,172 +1,172 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkEnergyComputer.h" #include #include using namespace mitk; EnergyComputer::EnergyComputer(ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen) : m_UseTrilinearInterpolation(true) { m_ParticleGrid = particleGrid; m_RandGen = randGen; //m_Image = qballImage; m_SphereInterpolator = interpolator; m_Mask = mask; m_ParticleLength = m_ParticleGrid->m_ParticleLength; m_SquaredParticleLength = m_ParticleLength*m_ParticleLength; m_Size[0] = mask->GetLargestPossibleRegion().GetSize()[0]; m_Size[1] = mask->GetLargestPossibleRegion().GetSize()[1]; m_Size[2] = mask->GetLargestPossibleRegion().GetSize()[2]; if (m_Size[0]<3 || m_Size[1]<3 || m_Size[2]<3) m_UseTrilinearInterpolation = false; m_Spacing[0] = mask->GetSpacing()[0]; m_Spacing[1] = mask->GetSpacing()[1]; m_Spacing[2] = mask->GetSpacing()[2]; // calculate rotation matrix vnl_matrix temp = mask->GetDirection().GetVnlMatrix(); vnl_matrix directionMatrix; directionMatrix.set_size(3,3); vnl_copy(temp, directionMatrix); vnl_vector_fixed d0 = directionMatrix.get_column(0); d0.normalize(); vnl_vector_fixed d1 = directionMatrix.get_column(1); d1.normalize(); vnl_vector_fixed d2 = directionMatrix.get_column(2); d2.normalize(); directionMatrix.set_column(0, d0); directionMatrix.set_column(1, d1); directionMatrix.set_column(2, d2); vnl_matrix_fixed I = directionMatrix*directionMatrix.transpose(); if(!I.is_identity(mitk::eps)) - fprintf(stderr,"itkGibbsTrackingFilter: image direction is not a rotation matrix. Tracking not possible!\n"); + fprintf(stderr,"EnergyComputer: image direction is not a rotation matrix. Tracking not possible!\n"); m_RotationMatrix = directionMatrix; if (QBALL_ODFSIZE != m_SphereInterpolator->nverts) fprintf(stderr,"EnergyComputer: error during init: data does not match with interpolation scheme\n"); int totsz = m_Size[0]*m_Size[1]*m_Size[2]; m_CumulatedSpatialProbability.resize(totsz + 1, 0.0); m_ActiveIndices.resize(totsz, 0); // calculate active voxels and cumulate probabilities m_NumActiveVoxels = 0; m_CumulatedSpatialProbability[0] = 0; for (int x = 0; x < m_Size[0];x++) for (int y = 0; y < m_Size[1];y++) for (int z = 0; z < m_Size[2];z++) { int idx = x+(y+z*m_Size[1])*m_Size[0]; ItkFloatImageType::IndexType index; index[0] = x; index[1] = y; index[2] = z; if (m_Mask->GetPixel(index) > 0.5) { m_CumulatedSpatialProbability[m_NumActiveVoxels+1] = m_CumulatedSpatialProbability[m_NumActiveVoxels] + m_Mask->GetPixel(index); m_ActiveIndices[m_NumActiveVoxels] = idx; m_NumActiveVoxels++; } } for (int k = 0; k < m_NumActiveVoxels; k++) m_CumulatedSpatialProbability[k] /= m_CumulatedSpatialProbability[m_NumActiveVoxels]; std::cout << "EnergyComputer: " << m_NumActiveVoxels << " active voxels found" << std::endl; } void EnergyComputer::SetParameters(float particleWeight, float particleWidth, float connectionPotential, float curvThres, float inexBalance, float particlePotential) { m_ParticleChemicalPotential = particlePotential; m_ConnectionPotential = connectionPotential; m_ParticleWeight = particleWeight; float bal = 1/(1+exp(-inexBalance)); m_ExtStrength = 2*bal; m_IntStrength = 2*(1-bal)/m_SquaredParticleLength; m_CurvatureThreshold = curvThres; float sigma_s = particleWidth; gamma_s = 1/(sigma_s*sigma_s); gamma_reg_s =1/(m_SquaredParticleLength/4); } // draw random position from active voxels void EnergyComputer::DrawRandomPosition(vnl_vector_fixed& R) { float r = m_RandGen->GetVariate();//m_RandGen->frand(); int j; int rl = 1; int rh = m_NumActiveVoxels; while(rh != rl) { j = rl + (rh-rl)/2; if (r < m_CumulatedSpatialProbability[j]) { rh = j; continue; } if (r > m_CumulatedSpatialProbability[j]) { rl = j+1; continue; } break; } R[0] = m_Spacing[0]*((float)(m_ActiveIndices[rh-1] % m_Size[0]) + m_RandGen->GetVariate()); R[1] = m_Spacing[1]*((float)((m_ActiveIndices[rh-1]/m_Size[0]) % m_Size[1]) + m_RandGen->GetVariate()); R[2] = m_Spacing[2]*((float)(m_ActiveIndices[rh-1]/(m_Size[0]*m_Size[1])) + m_RandGen->GetVariate()); } // return spatial probability of position float EnergyComputer::SpatProb(vnl_vector_fixed pos) { ItkFloatImageType::IndexType index; index[0] = floor(pos[0]/m_Spacing[0]); index[1] = floor(pos[1]/m_Spacing[1]); index[2] = floor(pos[2]/m_Spacing[2]); if (m_Mask->GetLargestPossibleRegion().IsInside(index)) // is inside image? return m_Mask->GetPixel(index); else return 0; } float EnergyComputer::mbesseli0(float x) { // BESSEL_APPROXCOEFF[0] = -0.1714; // BESSEL_APPROXCOEFF[1] = 0.5332; // BESSEL_APPROXCOEFF[2] = -1.4889; // BESSEL_APPROXCOEFF[3] = 2.0389; float y = x*x; float erg = -0.1714; erg += y*0.5332; erg += y*y*-1.4889; erg += y*y*y*2.0389; return erg; } float EnergyComputer::mexp(float x) { return((x>=7.0) ? 0 : ((x>=5.0) ? (-0.0029*x+0.0213) : ((x>=3.0) ? (-0.0215*x+0.1144) : ((x>=2.0) ? (-0.0855*x+0.3064) : ((x>=1.0) ? (-0.2325*x+0.6004) : ((x>=0.5) ? (-0.4773*x+0.8452) : ((x>=0.0) ? (-0.7869*x+1.0000) : 1 ))))))); // return exp(-x); } int EnergyComputer::GetNumActiveVoxels() { return m_NumActiveVoxels; } diff --git a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp index e6e4124b70..28f51c8bd2 100644 --- a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.cpp @@ -1,652 +1,699 @@ /*=================================================================== 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 __itkStreamlineTrackingFilter_txx #define __itkStreamlineTrackingFilter_txx #include #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_AngularThreshold(0.7), m_F(1.0), m_G(0.0), - m_Interpolate(true) + 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]); } std::cout << "StreamlineTrackingFilter: Angular threshold: " << m_AngularThreshold << 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; + if (m_Interpolate) + std::cout << "StreamlineTrackingFilter: using trilinear interpolation" << std::endl; + else + std::cout << "StreamlineTrackingFilter: using nearest neighbor interpolation" << 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> -void StreamlineTrackingFilter< TTensorPixelType, 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)) + 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; // 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()GetLargestPossibleRegion().IsInside(index) || m_FaImage->GetPixel(index)GetPixel(index)==0) + // is new position valid (inside image, above FA threshold etc.) + if (!IsValidPosition(pos, index, interpWeights)) // if not add last point and end streamline { m_InputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); - break; + return tractLength; + } + else if (distance>=m_PointPistance) + { + m_InputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); + ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); + distance = 0; } - if (!m_Interpolate) + 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 = dot_product(dirOld, dir); if (angleTransformContinuousIndexToPhysicalPoint( pos, worldPos ); - ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); } else dir = dirOld; } - else + else // use trilinear interpolation (weights calculated in IsValidPosition()) { - 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) - continue; - if (index[1] < 0 || index[1] >= m_ImageSize[1]-1) - continue; - if (index[2] < 0 || index[2] >= m_ImageSize[2]-1) - continue; - - typename InputImageType::IndexType tmpIdx; - typename InputImageType::PixelType tensor = m_InputImage->GetPixel(index); - tensor *= ( frac_x)*( frac_y)*( frac_z); - - tmpIdx = index; tmpIdx[0]++; - tensor += m_InputImage->GetPixel(tmpIdx) * (1-frac_x)*( frac_y)*( frac_z); + typename InputImageType::PixelType tensor = m_InputImage->GetPixel(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) * ( frac_x)*(1-frac_y)*( frac_z); + tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[2]; tmpIdx = index; tmpIdx[2]++; - tensor += m_InputImage->GetPixel(tmpIdx) * ( frac_x)*( frac_y)*(1-frac_z); + tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[3]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; - tensor += m_InputImage->GetPixel(tmpIdx) * (1-frac_x)*(1-frac_y)*( frac_z); + tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[4]; tmpIdx = index; tmpIdx[1]++; tmpIdx[2]++; - tensor += m_InputImage->GetPixel(tmpIdx) * ( frac_x)*(1-frac_y)*(1-frac_z); + tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[5]; tmpIdx = index; tmpIdx[2]++; tmpIdx[0]++; - tensor += m_InputImage->GetPixel(tmpIdx) * (1-frac_x)*( frac_y)*(1-frac_z); + tensor += m_InputImage->GetPixel(tmpIdx) * interpWeights[6]; tmpIdx = index; tmpIdx[0]++; tmpIdx[1]++; tmpIdx[2]++; - tensor += m_InputImage->GetPixel(tmpIdx) * (1-frac_x)*(1-frac_y)*(1-frac_z); + 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 = dot_product(dirOld, dir); if (angle=m_PointPistance) - { - m_InputImage->TransformContinuousIndexToPhysicalPoint( pos, worldPos ); - ids.push_back(points->InsertNextPoint(worldPos.GetDataPointer())); - distance = 0; - } } } + return tractLength; } template< class TTensorPixelType, class TPDPixelType> 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 - FollowStreamline(start, 1, points, pointIDs); + 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 - FollowStreamline(start, -1, points, pointIDs); + tractLength += FollowStreamline(start, -1, points, pointIDs); - // add ids to line counter += pointIDs.size(); + + if (tractLengthGetPointIds()->InsertNextId(pointIDs.at(i)); - if (counter>1) - Cells->InsertNextCell(line); + 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/Tractography/itkStreamlineTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h index 18b0bb2438..acbf55ee12 100644 --- a/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkStreamlineTrackingFilter.h @@ -1,133 +1,136 @@ /*=================================================================== 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. ===================================================================*/ /*=================================================================== This file is based heavily on a corresponding ITK filter. ===================================================================*/ #ifndef __itkStreamlineTrackingFilter_h_ #define __itkStreamlineTrackingFilter_h_ #include "MitkDiffusionImagingExports.h" #include #include #include #include #include #include #include #include #include namespace itk{ /** \class StreamlineTrackingFilter */ template< class TTensorPixelType, class TPDPixelType=double> class StreamlineTrackingFilter : public ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< Vector< TPDPixelType, 3 >, 3 > > { public: typedef StreamlineTrackingFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< Image< DiffusionTensor3D, 3 >, Image< Vector< TPDPixelType, 3 >, 3 > > Superclass; /** Method for creation through the object factory. */ itkNewMacro(Self) /** Runtime information support. */ itkTypeMacro(StreamlineTrackingFilter, ImageToImageFilter) typedef TTensorPixelType TensorComponentType; typedef TPDPixelType DirectionPixelType; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef itk::Image ItkUcharImgType; typedef itk::Image ItkFloatImgType; typedef itk::Image< vnl_vector_fixed, 3> ItkPDImgType; typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; itkGetMacro( FiberPolyData, FiberPolyDataType ) itkSetMacro( SeedImage, ItkUcharImgType::Pointer) itkSetMacro( MaskImage, ItkUcharImgType::Pointer) itkSetMacro( SeedsPerVoxel, int) itkSetMacro( FaThreshold, float) itkSetMacro( AngularThreshold, float) itkSetMacro( StepSize, float) itkSetMacro( F, float ) itkSetMacro( G, float ) itkSetMacro( Interpolate, bool ) + itkSetMacro( MinTractLength, float ) + itkGetMacro( MinTractLength, float ) protected: StreamlineTrackingFilter(); ~StreamlineTrackingFilter() {} void PrintSelf(std::ostream& os, Indent indent) const; void CalculateNewPosition(itk::ContinuousIndex& pos, vnl_vector_fixed& dir, typename InputImageType::IndexType& index); - void FollowStreamline(itk::ContinuousIndex pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids); - + float FollowStreamline(itk::ContinuousIndex pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids); + bool IsValidPosition(itk::ContinuousIndex& pos, typename InputImageType::IndexType& index, vnl_vector_fixed< float, 8 >& interpWeights); double RoundToNearest(double num); void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, int threadId); void AfterThreadedGenerateData(); FiberPolyDataType AddPolyData(FiberPolyDataType poly1, FiberPolyDataType poly2); FiberPolyDataType m_FiberPolyData; vtkSmartPointer m_Points; vtkSmartPointer m_Cells; ItkFloatImgType::Pointer m_EmaxImage; ItkFloatImgType::Pointer m_FaImage; ItkPDImgType::Pointer m_PdImage; typename InputImageType::Pointer m_InputImage; float m_FaThreshold; float m_AngularThreshold; float m_StepSize; int m_MaxLength; + float m_MinTractLength; int m_SeedsPerVoxel; float m_F; float m_G; std::vector< int > m_ImageSize; std::vector< float > m_ImageSpacing; ItkUcharImgType::Pointer m_SeedImage; ItkUcharImgType::Pointer m_MaskImage; bool m_Interpolate; float m_PointPistance; itk::VectorContainer< int, FiberPolyDataType >::Pointer m_PolyDataContainer; private: }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkStreamlineTrackingFilter.cpp" #endif #endif //__itkStreamlineTrackingFilter_h_ diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui index 17481a4b61..94dd852334 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkFiberProcessingViewControls.ui @@ -1,912 +1,912 @@ QmitkFiberProcessingViewControls 0 0 665 683 Form 0 9 3 9 3 Fiber Extraction 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 30 30 Draw circular ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/circle.png:/QmitkDiffusionImaging/circle.png 32 32 false true 30 30 Draw rectangular ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/rectangle.png:/QmitkDiffusionImaging/rectangle.png 32 32 true true 30 30 Draw polygonal ROI. Select reference fiber bundle to execute. :/QmitkDiffusionImaging/polygon.png:/QmitkDiffusionImaging/polygon.png 32 32 true true Qt::Horizontal 40 20 QFrame::NoFrame QFrame::Raised 0 false 0 0 200 16777215 11 Extract fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute. Extract false 0 0 200 16777215 11 Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute. Substract false 0 0 200 16777215 11 Merge selected fiber bundles. Select at least two fiber bundles to execute. Join Qt::Horizontal 40 20 false 0 0 200 16777215 11 Extract fibers passing through selected surface mesh. Select surface mesh and fiber bundle to execute. Extract 3D false 0 0 16777215 16777215 11 Generate a binary image containing all selected ROIs. Select at least one ROI (planar figure) and a reference fiber bundle or image. ROI Image 0 0 200 0 16777215 60 QFrame::NoFrame QFrame::Raised 0 Qt::Horizontal 40 20 false 60 16777215 Create AND composition with selected ROIs. AND false 60 16777215 Create OR composition with selected ROIs. OR false 60 16777215 Create NOT composition from selected ROI. NOT Fiber Processing QFormLayout::AllNonFixedFieldsGrow 0 0 Tract Density Image Tract Density Image (normalize image values) Binary Envelope Fiber Bundle Image Fiber Endings Image Fiber Endings Pointset false 0 0 200 16777215 11 Perform selected operation on all selected fiber bundles. Generate If selected operation generates an image, the inverse image is returned. Invert false 0 0 200 16777215 11 Resample fibers using a Kochanek spline interpolation. Smooth Fibers Points per cm 1 50 5 false 0 0 200 16777215 11 Remove fibers shorter/longer than the specified length (in mm). Length Threshold false 0 0 200 16777215 11 Mirror fibers around specified axis. Mirror Fibers false 0 0 200 16777215 11 Apply float image values (0-1) as color coding to the selected fiber bundle. Color By Scalar Map 0 3 3 Sagittal Coronal Transversal Upsampling factor 1 0.100000000000000 10.000000000000000 0.100000000000000 1.000000000000000 QFrame::NoFrame QFrame::Raised 0 0 Minimum fiber length in mm 0 1000 20 Maximum fiber length in mm 0 - 1000 + 10000 500 false 0 0 200 16777215 11 Remove fibers with a too high curvature Curvature Threshold QFrame::NoFrame QFrame::Raised 0 0 Minimum radius of circle created by three consecutive points of a fiber 100.000000000000000 0.100000000000000 2.000000000000000 Remove whole fiber if it is exceeding the curvature threshold, otherwise remove only high curvature part. Remove Fiber true Fiber Statistics Courier 10 Pitch false true Qt::Vertical 20 40 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui index 4b2ea27ef5..5bcbf87bbc 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkQBallReconstructionViewControls.ui @@ -1,282 +1,282 @@ QmitkQBallReconstructionViewControls 0 0 350 844 0 0 true QmitkQBallReconstructionViewControls Data Diffusion Image: - Reconstruction Advanced Settings QFrame::StyledPanel QFrame::Raised QFormLayout::AllNonFixedFieldsGrow true B0 Threshold false true 0 true Output B0-Image true Spherical Harmonics: true Maximum l-Level false true -1 true Regularization Parameter Lambda false true 0.006 Output SH-Coeffs-Image - 0 + 2 Numerical Standard (SH) Solid Angle (SH) Constraint Solid Angle (SH) ADC-Profile only (SH) Raw Signal only (SH) - Mulit q-Ball (SH) + Multi-Shell (SH) TextLabel false Start Reconstruction true Qt::LeftToRight false - Multi q-Ball reconstruction + Multi-Shell Reconstruction Qt::Vertical 20 0 diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.cpp index 8dd2b1e762..4b283b73c2 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingView.cpp @@ -1,241 +1,239 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include #include // Qmitk #include "QmitkStreamlineTrackingView.h" #include "QmitkStdMultiWidget.h" // Qt #include // MITK #include #include #include // VTK #include #include #include #include #include #include const std::string QmitkStreamlineTrackingView::VIEW_ID = "org.mitk.views.streamlinetracking"; const std::string id_DataManager = "org.mitk.views.datamanager"; using namespace berry; QmitkStreamlineTrackingView::QmitkStreamlineTrackingView() : QmitkFunctionality() , m_Controls( 0 ) , m_MultiWidget( NULL ) , m_TensorImage( NULL ) , m_SeedRoi( NULL ) { } // Destructor QmitkStreamlineTrackingView::~QmitkStreamlineTrackingView() { } void QmitkStreamlineTrackingView::CreateQtPartControl( QWidget *parent ) { if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkStreamlineTrackingViewControls; m_Controls->setupUi( parent ); connect( m_Controls->commandLinkButton, SIGNAL(clicked()), this, SLOT(DoFiberTracking()) ); connect( m_Controls->m_SeedsPerVoxelSlider, SIGNAL(valueChanged(int)), this, SLOT(OnSeedsPerVoxelChanged(int)) ); connect( m_Controls->m_MinTractLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(OnMinTractLengthChanged(int)) ); connect( m_Controls->m_FaThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(OnFaThresholdChanged(int)) ); connect( m_Controls->m_AngularThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(OnAngularThresholdChanged(int)) ); connect( m_Controls->m_StepsizeSlider, SIGNAL(valueChanged(int)), this, SLOT(OnStepsizeChanged(int)) ); connect( m_Controls->m_fSlider, SIGNAL(valueChanged(int)), this, SLOT(OnfChanged(int)) ); connect( m_Controls->m_gSlider, SIGNAL(valueChanged(int)), this, SLOT(OngChanged(int)) ); } } void QmitkStreamlineTrackingView::OnfChanged(int value) { m_Controls->m_fLabel->setText(QString("f: ")+QString::number((float)value/100)); } void QmitkStreamlineTrackingView::OngChanged(int value) { m_Controls->m_gLabel->setText(QString("g: ")+QString::number((float)value/100)); } void QmitkStreamlineTrackingView::OnAngularThresholdChanged(int value) { m_Controls->m_AngularThresholdLabel->setText(QString("Angular Threshold: ")+QString::number(value)+QString("°")); } void QmitkStreamlineTrackingView::OnSeedsPerVoxelChanged(int value) { m_Controls->m_SeedsPerVoxelLabel->setText(QString("Seeds per Voxel: ")+QString::number(value)); } void QmitkStreamlineTrackingView::OnMinTractLengthChanged(int value) { m_Controls->m_MinTractLengthLabel->setText(QString("Min. Tract Length: ")+QString::number(value)+QString("mm")); } void QmitkStreamlineTrackingView::OnFaThresholdChanged(int value) { m_Controls->m_FaThresholdLabel->setText(QString("FA Threshold: ")+QString::number((float)value/100)); } void QmitkStreamlineTrackingView::OnStepsizeChanged(int value) { if (value==0) m_Controls->m_StepsizeLabel->setText(QString("Stepsize: auto")); else m_Controls->m_StepsizeLabel->setText(QString("Stepsize: ")+QString::number((float)value/10)+QString("mm")); } void QmitkStreamlineTrackingView::StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkStreamlineTrackingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } void QmitkStreamlineTrackingView::OnSelectionChanged( std::vector nodes ) { m_TensorImageNode = NULL; m_TensorImage = NULL; m_SeedRoi = NULL; m_MaskImage = NULL; m_Controls->m_TensorImageLabel->setText("-"); m_Controls->m_RoiImageLabel->setText("-"); m_Controls->m_MaskImageLabel->setText("-"); if(nodes.empty()) return; for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { if( dynamic_cast(node->GetData()) ) { m_TensorImageNode = node; m_TensorImage = dynamic_cast(node->GetData()); m_Controls->m_TensorImageLabel->setText(node->GetName().c_str()); } else { bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary && m_SeedRoi.IsNull()) { m_SeedRoi = dynamic_cast(node->GetData()); m_Controls->m_RoiImageLabel->setText(node->GetName().c_str()); } else if (isBinary) { m_MaskImage = dynamic_cast(node->GetData()); m_Controls->m_MaskImageLabel->setText(node->GetName().c_str()); } } } } if(m_TensorImage.IsNotNull()) m_Controls->commandLinkButton->setEnabled(true); else m_Controls->commandLinkButton->setEnabled(false); } void QmitkStreamlineTrackingView::DoFiberTracking() { if (m_TensorImage.IsNull()) return; typedef itk::Image< itk::DiffusionTensor3D, 3> TensorImageType; typedef mitk::ImageToItk CastType; typedef mitk::ImageToItk CastType2; CastType::Pointer caster = CastType::New(); caster->SetInput(m_TensorImage); caster->Update(); TensorImageType::Pointer image = caster->GetOutput(); typedef itk::StreamlineTrackingFilter< float > FilterType; FilterType::Pointer filter = FilterType::New(); filter->SetInput(image); filter->SetSeedsPerVoxel(m_Controls->m_SeedsPerVoxelSlider->value()); filter->SetFaThreshold((float)m_Controls->m_FaThresholdSlider->value()/100); filter->SetAngularThreshold(cos((float)m_Controls->m_AngularThresholdSlider->value()*M_PI/180)); filter->SetStepSize((float)m_Controls->m_StepsizeSlider->value()/10); filter->SetF((float)m_Controls->m_fSlider->value()/100); filter->SetG((float)m_Controls->m_gSlider->value()/100); filter->SetInterpolate(m_Controls->m_InterpolationBox->isChecked()); + filter->SetMinTractLength(m_Controls->m_MinTractLengthSlider->value()); //filter->SetNumberOfThreads(1); if (m_SeedRoi.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(m_SeedRoi, mask); filter->SetSeedImage(mask); } if (m_MaskImage.IsNotNull()) { ItkUCharImageType::Pointer mask = ItkUCharImageType::New(); mitk::CastToItkImage(m_MaskImage, mask); filter->SetMaskImage(mask); } filter->Update(); vtkSmartPointer fiberBundle = filter->GetFiberPolyData(); if ( fiberBundle->GetNumberOfLines()==0 ) return; mitk::FiberBundleX::Pointer fib = mitk::FiberBundleX::New(fiberBundle); - if (fib->RemoveShortFibers(m_Controls->m_MinTractLengthSlider->value())) - { - mitk::DataNode::Pointer node = mitk::DataNode::New(); - node->SetData(fib); - QString name(m_TensorImageNode->GetName().c_str()); - name += "_FiberBundle"; - node->SetName(name.toStdString()); - node->SetVisibility(true); - GetDataStorage()->Add(node); - } + mitk::DataNode::Pointer node = mitk::DataNode::New(); + node->SetData(fib); + QString name(m_TensorImageNode->GetName().c_str()); + name += "_Streamline"; + node->SetName(name.toStdString()); + node->SetVisibility(true); + GetDataStorage()->Add(node); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui index 64c86d2904..aa5e2fb2bf 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkStreamlineTrackingViewControls.ui @@ -1,385 +1,385 @@ QmitkStreamlineTrackingViewControls 0 0 480 553 0 0 QmitkTemplate 3 3 0 Qt::Vertical QSizePolicy::Expanding 20 220 Parameters Number of tracts started in each voxel of the seed ROI. 1 - 10 + 100 Qt::Horizontal Fractional Anisotropy Threshold 0 90 45 Qt::Horizontal Minimum tract length in mm. 0 500 40 Qt::Horizontal Number of tracts started in each voxel of the seed ROI. Seeds per Voxel: 1 Qt::Horizontal QSizePolicy::Fixed 200 0 Fractional Anisotropy Threshold 0 100 20 Qt::Horizontal Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0). 0 100 100 Qt::Horizontal Minimum tract length in mm. FA Threshold: 0.2 Minimum tract length in mm. Min. Tract Length: 40mm Stepsize in mm (auto = 0.1*minimal spacing) 0 100 0 Qt::Horizontal Minimum tract length in mm. g: 0 Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking) 0 100 0 Qt::Horizontal Minimum tract length in mm. Angular Threshold: 45° Minimum tract length in mm. Step Size: auto Minimum tract length in mm. f: 1 Default is nearest neighbor interpolation. Enable trilinear interpolation Number of tracts started in each voxel of the seed ROI. Mori et al. Annals Neurology 1999 false Start Tracking Data - - Seed ROI Image: Tensor Image: Stop tracking if leaving mask Mask Image: - Number of tracts started in each voxel of the seed ROI. Lazar et al. Human Brain Mapping 2003 Number of tracts started in each voxel of the seed ROI. Weinstein et al. Proceedings of IEEE Visualization 1999 commandLinkButton