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 <vnl/vnl_copy.h>
#include <itkNumericTraits.h>
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<double> temp = mask->GetDirection().GetVnlMatrix();
vnl_matrix<float> directionMatrix; directionMatrix.set_size(3,3);
vnl_copy(temp, directionMatrix);
vnl_vector_fixed<float, 3> d0 = directionMatrix.get_column(0); d0.normalize();
vnl_vector_fixed<float, 3> d1 = directionMatrix.get_column(1); d1.normalize();
vnl_vector_fixed<float, 3> 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<float, 3, 3> 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<float, 3>& 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<float, 3> 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 <time.h>
#include <stdio.h>
#include <stdlib.h>
#include "itkStreamlineTrackingFilter.h"
#include <itkImageRegionConstIterator.h>
#include <itkImageRegionConstIteratorWithIndex.h>
#include <itkImageRegionIterator.h>
#define _USE_MATH_DEFINES
#include <math.h>
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]<m_ImageSpacing[1] && m_ImageSpacing[0]<m_ImageSpacing[2])
minSpacing = m_ImageSpacing[0];
else if (m_ImageSpacing[1] < m_ImageSpacing[2])
minSpacing = m_ImageSpacing[1];
else
minSpacing = m_ImageSpacing[2];
if (m_StepSize<0.1*minSpacing)
{
m_StepSize = 0.1*minSpacing;
m_PointPistance = 0.5*minSpacing;
}
m_PolyDataContainer = itk::VectorContainer< int, FiberPolyDataType >::New();
for (int i=0; i<this->GetNumberOfThreads(); 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<TTensorPixelType> TensorType;
typename TensorType::EigenValuesArrayType eigenvalues;
typename TensorType::EigenVectorsMatrixType eigenvectors;
for (int x=0; x<m_ImageSize[0]; x++)
for (int y=0; y<m_ImageSize[1]; y++)
for (int z=0; z<m_ImageSize[2]; z++)
{
typename InputImageType::IndexType index;
index[0] = x; index[1] = y; index[2] = z;
typename InputImageType::PixelType tensor = m_InputImage->GetPixel(index);
vnl_vector_fixed<double,3> 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<double, 3>& pos, vnl_vector_fixed<double,3>& 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<double, 3>& 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 (FA<m_FaThreshold)
+ return false;
+ }
+ else if (m_FaImage->GetPixel(index)<m_FaThreshold)
+ return false;
+
+ return true;
+}
+
+template< class TTensorPixelType, class TPDPixelType>
+float StreamlineTrackingFilter< TTensorPixelType, TPDPixelType>
::FollowStreamline(itk::ContinuousIndex<double, 3> pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids)
{
+ float tractLength = 0;
typedef itk::DiffusionTensor3D<TTensorPixelType> 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<double> 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<double,3> dir = m_PdImage->GetPixel(index);
dir *= dirSign; // reverse direction
vnl_vector_fixed<double,3> dirOld = dir;
if (dir.magnitude()<mitk::eps)
- return;
+ return tractLength;
for (int step=0; step< m_MaxLength/2; step++)
{
// get new position
CalculateNewPosition(pos, dir, index);
distance += m_StepSize;
+ tractLength += m_StepSize;
- // are we still inside the image and is the fa value high enough?
- if (!m_InputImage->GetLargestPossibleRegion().IsInside(index) || m_FaImage->GetPixel(index)<m_FaThreshold || m_MaskImage->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 (angle<m_AngularThreshold)
break;
if (dir.magnitude()<mitk::eps)
dir = dirOld;
else
dirOld = dir;
indexOld = index;
-
- // add position to streamline
- m_InputImage->TransformContinuousIndexToPhysicalPoint( 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_AngularThreshold)
break;
if (dir.magnitude()<mitk::eps)
dir = dirOld;
else
dirOld = dir;
indexOld = index;
-
- // add position to streamline
- if (distance>=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<vtkPoints> points = vtkPoints::New();
vtkSmartPointer<vtkCellArray> Cells = vtkCellArray::New();
typedef itk::DiffusionTensor3D<TTensorPixelType> 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<double> worldPos;
while( !it.IsAtEnd() )
{
if (mit.Value()==0 || fit.Value()<m_FaThreshold || mit2.Value()==0)
{
++mit;
++mit2;
++it;
++fit;
continue;
}
for (int s=0; s<m_SeedsPerVoxel; s++)
{
vtkSmartPointer<vtkPolyLine> line = vtkSmartPointer<vtkPolyLine>::New();
std::vector< vtkIdType > pointIDs;
typename InputImageType::IndexType index = it.GetIndex();
itk::ContinuousIndex<double, 3> 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 (tractLength<m_MinTractLength || counter<2)
+ continue;
+
+ // add ids to line
for (int i=0; i<pointIDs.size(); i++)
line->GetPointIds()->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<vtkPolyData> vNewPolyData = vtkSmartPointer<vtkPolyData>::New();
vtkSmartPointer<vtkCellArray> vNewLines = poly1->GetLines();
vtkSmartPointer<vtkPoints> vNewPoints = poly1->GetPoints();
vtkSmartPointer<vtkCellArray> vLines = poly2->GetLines();
vLines->InitTraversal();
for( int i=0; i<vLines->GetNumberOfCells(); i++ )
{
vtkIdType numPoints(0);
vtkIdType* points(NULL);
vLines->GetNextCell ( numPoints, points );
vtkSmartPointer<vtkPolyLine> container = vtkSmartPointer<vtkPolyLine>::New();
for( int j=0; j<numPoints; j++)
{
vtkIdType id = vNewPoints->InsertNextPoint(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; i<this->GetNumberOfThreads(); 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 <itkImageToImageFilter.h>
#include <itkVectorContainer.h>
#include <itkVectorImage.h>
#include <itkDiffusionTensor3D.h>
#include <vtkSmartPointer.h>
#include <vtkPolyData.h>
#include <vtkCellArray.h>
#include <vtkPoints.h>
#include <vtkPolyLine.h>
namespace itk{
/** \class StreamlineTrackingFilter
*/
template< class TTensorPixelType, class TPDPixelType=double>
class StreamlineTrackingFilter :
public ImageToImageFilter< Image< DiffusionTensor3D<TTensorPixelType>, 3 >,
Image< Vector< TPDPixelType, 3 >, 3 > >
{
public:
typedef StreamlineTrackingFilter Self;
typedef SmartPointer<Self> Pointer;
typedef SmartPointer<const Self> ConstPointer;
typedef ImageToImageFilter< Image< DiffusionTensor3D<TTensorPixelType>, 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<unsigned char, 3> ItkUcharImgType;
typedef itk::Image<float, 3> ItkFloatImgType;
typedef itk::Image< vnl_vector_fixed<double,3>, 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<double, 3>& pos, vnl_vector_fixed<double,3>& dir, typename InputImageType::IndexType& index);
- void FollowStreamline(itk::ContinuousIndex<double, 3> pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids);
-
+ float FollowStreamline(itk::ContinuousIndex<double, 3> pos, int dirSign, vtkPoints* points, std::vector< vtkIdType >& ids);
+ bool IsValidPosition(itk::ContinuousIndex<double, 3>& 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<vtkPoints> m_Points;
vtkSmartPointer<vtkCellArray> 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 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>QmitkFiberProcessingViewControls</class>
<widget class="QWidget" name="QmitkFiberProcessingViewControls">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>665</width>
<height>683</height>
</rect>
</property>
<property name="windowTitle">
<string>Form</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout">
<property name="spacing">
<number>0</number>
</property>
<property name="leftMargin">
<number>9</number>
</property>
<property name="topMargin">
<number>3</number>
</property>
<property name="rightMargin">
<number>9</number>
</property>
<property name="bottomMargin">
<number>3</number>
</property>
<item>
<widget class="QGroupBox" name="groupBox">
<property name="title">
<string>Fiber Extraction</string>
</property>
<layout class="QGridLayout" name="gridLayout_4">
<item row="0" column="0">
<widget class="QFrame" name="m_PlanarFigureButtonsFrame">
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="minimumSize">
<size>
<width>200</width>
<height>0</height>
</size>
</property>
<property name="maximumSize">
<size>
<width>16777215</width>
<height>60</height>
</size>
</property>
<property name="frameShape">
<enum>QFrame::NoFrame</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
<layout class="QGridLayout" name="gridLayout">
<property name="margin">
<number>0</number>
</property>
<item row="0" column="0">
<widget class="QPushButton" name="m_CircleButton">
<property name="maximumSize">
<size>
<width>30</width>
<height>30</height>
</size>
</property>
<property name="toolTip">
<string>Draw circular ROI. Select reference fiber bundle to execute.</string>
</property>
<property name="text">
<string/>
</property>
<property name="icon">
<iconset resource="../../resources/QmitkDiffusionImaging.qrc">
<normaloff>:/QmitkDiffusionImaging/circle.png</normaloff>:/QmitkDiffusionImaging/circle.png</iconset>
</property>
<property name="iconSize">
<size>
<width>32</width>
<height>32</height>
</size>
</property>
<property name="checkable">
<bool>false</bool>
</property>
<property name="flat">
<bool>true</bool>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QPushButton" name="m_RectangleButton">
<property name="maximumSize">
<size>
<width>30</width>
<height>30</height>
</size>
</property>
<property name="toolTip">
<string>Draw rectangular ROI. Select reference fiber bundle to execute.</string>
</property>
<property name="text">
<string/>
</property>
<property name="icon">
<iconset resource="../../resources/QmitkDiffusionImaging.qrc">
<normaloff>:/QmitkDiffusionImaging/rectangle.png</normaloff>:/QmitkDiffusionImaging/rectangle.png</iconset>
</property>
<property name="iconSize">
<size>
<width>32</width>
<height>32</height>
</size>
</property>
<property name="checkable">
<bool>true</bool>
</property>
<property name="flat">
<bool>true</bool>
</property>
</widget>
</item>
<item row="0" column="2">
<widget class="QPushButton" name="m_PolygonButton">
<property name="maximumSize">
<size>
<width>30</width>
<height>30</height>
</size>
</property>
<property name="toolTip">
<string>Draw polygonal ROI. Select reference fiber bundle to execute.</string>
</property>
<property name="text">
<string/>
</property>
<property name="icon">
<iconset resource="../../resources/QmitkDiffusionImaging.qrc">
<normaloff>:/QmitkDiffusionImaging/polygon.png</normaloff>:/QmitkDiffusionImaging/polygon.png</iconset>
</property>
<property name="iconSize">
<size>
<width>32</width>
<height>32</height>
</size>
</property>
<property name="checkable">
<bool>true</bool>
</property>
<property name="flat">
<bool>true</bool>
</property>
</widget>
</item>
<item row="0" column="3">
<spacer name="horizontalSpacer">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>40</width>
<height>20</height>
</size>
</property>
</spacer>
</item>
</layout>
</widget>
</item>
<item row="2" column="0">
<widget class="QFrame" name="frame_2">
<property name="frameShape">
<enum>QFrame::NoFrame</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
<layout class="QGridLayout" name="gridLayout_3">
<property name="margin">
<number>0</number>
</property>
<item row="0" column="0">
<widget class="QCommandLinkButton" name="doExtractFibersButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Extract fibers passing through selected ROI or composite ROI. Select ROI and fiber bundle to execute.</string>
</property>
<property name="text">
<string>Extract</string>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QCommandLinkButton" name="m_SubstractBundles">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Returns all fibers contained in bundle X that are not contained in bundle Y (not commutative!). Select at least two fiber bundles to execute.</string>
</property>
<property name="text">
<string>Substract</string>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="QCommandLinkButton" name="m_JoinBundles">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Merge selected fiber bundles. Select at least two fiber bundles to execute.</string>
</property>
<property name="text">
<string>Join</string>
</property>
</widget>
</item>
<item row="0" column="2">
<spacer name="horizontalSpacer_2">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>40</width>
<height>20</height>
</size>
</property>
</spacer>
</item>
<item row="0" column="1">
<widget class="QCommandLinkButton" name="m_Extract3dButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Extract fibers passing through selected surface mesh. Select surface mesh and fiber bundle to execute.</string>
</property>
<property name="text">
<string>Extract 3D</string>
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QCommandLinkButton" name="m_GenerateRoiImage">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>16777215</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Generate a binary image containing all selected ROIs. Select at least one ROI (planar figure) and a reference fiber bundle or image.</string>
</property>
<property name="text">
<string>ROI Image</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item row="1" column="0">
<widget class="QFrame" name="m_PlanarFigureButtonsFrame_2">
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="minimumSize">
<size>
<width>200</width>
<height>0</height>
</size>
</property>
<property name="maximumSize">
<size>
<width>16777215</width>
<height>60</height>
</size>
</property>
<property name="frameShape">
<enum>QFrame::NoFrame</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
<layout class="QGridLayout" name="gridLayout_5">
<property name="margin">
<number>0</number>
</property>
<item row="0" column="3">
<spacer name="horizontalSpacer_3">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>40</width>
<height>20</height>
</size>
</property>
</spacer>
</item>
<item row="0" column="0">
<widget class="QCommandLinkButton" name="PFCompoANDButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="maximumSize">
<size>
<width>60</width>
<height>16777215</height>
</size>
</property>
<property name="toolTip">
<string>Create AND composition with selected ROIs.</string>
</property>
<property name="text">
<string>AND</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QCommandLinkButton" name="PFCompoORButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="maximumSize">
<size>
<width>60</width>
<height>16777215</height>
</size>
</property>
<property name="toolTip">
<string>Create OR composition with selected ROIs.</string>
</property>
<property name="text">
<string>OR</string>
</property>
</widget>
</item>
<item row="0" column="2">
<widget class="QCommandLinkButton" name="PFCompoNOTButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="maximumSize">
<size>
<width>60</width>
<height>16777215</height>
</size>
</property>
<property name="toolTip">
<string>Create NOT composition from selected ROI.</string>
</property>
<property name="text">
<string>NOT</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<widget class="QGroupBox" name="groupBox_2">
<property name="title">
<string>Fiber Processing</string>
</property>
<layout class="QFormLayout" name="formLayout">
<property name="fieldGrowthPolicy">
<enum>QFormLayout::AllNonFixedFieldsGrow</enum>
</property>
<item row="0" column="0">
<widget class="QComboBox" name="m_GenerationBox">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Fixed">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<item>
<property name="text">
<string>Tract Density Image</string>
</property>
</item>
<item>
<property name="text">
<string>Tract Density Image (normalize image values)</string>
</property>
</item>
<item>
<property name="text">
<string>Binary Envelope</string>
</property>
</item>
<item>
<property name="text">
<string>Fiber Bundle Image</string>
</property>
</item>
<item>
<property name="text">
<string>Fiber Endings Image</string>
</property>
</item>
<item>
<property name="text">
<string>Fiber Endings Pointset</string>
</property>
</item>
</widget>
</item>
<item row="1" column="0">
<widget class="QCommandLinkButton" name="m_ProcessFiberBundleButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Perform selected operation on all selected fiber bundles.</string>
</property>
<property name="text">
<string>Generate</string>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="QCheckBox" name="m_InvertCheckbox">
<property name="toolTip">
<string>If selected operation generates an image, the inverse image is returned.</string>
</property>
<property name="text">
<string>Invert</string>
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QCommandLinkButton" name="m_ResampleFibersButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Resample fibers using a Kochanek spline interpolation.</string>
</property>
<property name="text">
<string>Smooth Fibers</string>
</property>
</widget>
</item>
<item row="2" column="1">
<widget class="QSpinBox" name="m_ResampleFibersSpinBox">
<property name="toolTip">
<string>Points per cm</string>
</property>
<property name="minimum">
<number>1</number>
</property>
<property name="maximum">
<number>50</number>
</property>
<property name="value">
<number>5</number>
</property>
</widget>
</item>
<item row="3" column="0">
<widget class="QCommandLinkButton" name="m_PruneFibersButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Remove fibers shorter/longer than the specified length (in mm).</string>
</property>
<property name="text">
<string>Length Threshold</string>
</property>
</widget>
</item>
<item row="6" column="0">
<widget class="QCommandLinkButton" name="m_MirrorFibersButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Mirror fibers around specified axis.</string>
</property>
<property name="text">
<string>Mirror Fibers</string>
</property>
</widget>
</item>
<item row="7" column="0">
<widget class="QCommandLinkButton" name="m_FaColorFibersButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Apply float image values (0-1) as color coding to the selected fiber bundle.</string>
</property>
<property name="text">
<string>Color By Scalar Map</string>
</property>
</widget>
</item>
<item row="6" column="1">
<widget class="QComboBox" name="m_AxisSelectionBox">
<property name="currentIndex">
<number>0</number>
</property>
<property name="maxVisibleItems">
<number>3</number>
</property>
<property name="maxCount">
<number>3</number>
</property>
<item>
<property name="text">
<string>Sagittal</string>
</property>
</item>
<item>
<property name="text">
<string>Coronal</string>
</property>
</item>
<item>
<property name="text">
<string>Transversal</string>
</property>
</item>
</widget>
</item>
<item row="0" column="1">
<widget class="QDoubleSpinBox" name="m_UpsamplingSpinBox">
<property name="toolTip">
<string>Upsampling factor</string>
</property>
<property name="decimals">
<number>1</number>
</property>
<property name="minimum">
<double>0.100000000000000</double>
</property>
<property name="maximum">
<double>10.000000000000000</double>
</property>
<property name="singleStep">
<double>0.100000000000000</double>
</property>
<property name="value">
<double>1.000000000000000</double>
</property>
</widget>
</item>
<item row="3" column="1">
<widget class="QFrame" name="frame">
<property name="frameShape">
<enum>QFrame::NoFrame</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
<layout class="QGridLayout" name="gridLayout_6">
<property name="margin">
<number>0</number>
</property>
<property name="spacing">
<number>0</number>
</property>
<item row="0" column="0">
<widget class="QSpinBox" name="m_PruneFibersSpinBox">
<property name="toolTip">
<string>Minimum fiber length in mm</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
<number>1000</number>
</property>
<property name="value">
<number>20</number>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QSpinBox" name="m_MaxPruneFibersSpinBox">
<property name="toolTip">
<string>Maximum fiber length in mm</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
- <number>1000</number>
+ <number>10000</number>
</property>
<property name="value">
<number>500</number>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item row="4" column="0">
<widget class="QCommandLinkButton" name="m_CurvatureThresholdButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="sizePolicy">
<sizepolicy hsizetype="Preferred" vsizetype="Preferred">
<horstretch>0</horstretch>
<verstretch>0</verstretch>
</sizepolicy>
</property>
<property name="maximumSize">
<size>
<width>200</width>
<height>16777215</height>
</size>
</property>
<property name="font">
<font>
<pointsize>11</pointsize>
</font>
</property>
<property name="toolTip">
<string>Remove fibers with a too high curvature</string>
</property>
<property name="text">
<string>Curvature Threshold</string>
</property>
</widget>
</item>
<item row="4" column="1">
<widget class="QFrame" name="frame_3">
<property name="frameShape">
<enum>QFrame::NoFrame</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
<layout class="QGridLayout" name="gridLayout_7">
<property name="margin">
<number>0</number>
</property>
<property name="spacing">
<number>0</number>
</property>
<item row="0" column="0">
<widget class="QDoubleSpinBox" name="m_MinCurvatureRadiusBox">
<property name="toolTip">
<string>Minimum radius of circle created by three consecutive points of a fiber</string>
</property>
<property name="maximum">
<double>100.000000000000000</double>
</property>
<property name="singleStep">
<double>0.100000000000000</double>
</property>
<property name="value">
<double>2.000000000000000</double>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QCheckBox" name="m_RemoveFiberDueToCurvatureCheckbox">
<property name="toolTip">
<string>Remove whole fiber if it is exceeding the curvature threshold, otherwise remove only high curvature part.</string>
</property>
<property name="text">
<string>Remove Fiber</string>
</property>
<property name="checked">
<bool>true</bool>
</property>
</widget>
</item>
</layout>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<widget class="QGroupBox" name="groupBox_3">
<property name="title">
<string>Fiber Statistics</string>
</property>
<layout class="QGridLayout" name="gridLayout_2">
<item row="0" column="0">
<widget class="QTextEdit" name="m_StatsTextEdit">
<property name="font">
<font>
<family>Courier 10 Pitch</family>
</font>
</property>
<property name="acceptDrops">
<bool>false</bool>
</property>
<property name="readOnly">
<bool>true</bool>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<spacer name="verticalSpacer">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
</property>
</spacer>
</item>
</layout>
</widget>
<resources>
<include location="../../resources/QmitkDiffusionImaging.qrc"/>
</resources>
<connections/>
</ui>
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 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>QmitkQBallReconstructionViewControls</class>
<widget class="QWidget" name="QmitkQBallReconstructionViewControls">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>350</width>
<height>844</height>
</rect>
</property>
<property name="minimumSize">
<size>
<width>0</width>
<height>0</height>
</size>
</property>
<property name="acceptDrops">
<bool>true</bool>
</property>
<property name="windowTitle">
<string>QmitkQBallReconstructionViewControls</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout_3">
<item>
<widget class="QGroupBox" name="groupBox">
<property name="title">
<string>Data</string>
</property>
<layout class="QGridLayout" name="gridLayout">
<item row="0" column="0">
<widget class="QLabel" name="label">
<property name="text">
<string>Diffusion Image:</string>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QLabel" name="m_DiffusionImageLabel">
<property name="text">
<string>-</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<widget class="QGroupBox" name="groupBox_3">
<property name="title">
<string>Reconstruction</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout">
<item>
<widget class="QCheckBox" name="m_AdvancedCheckbox">
<property name="text">
<string>Advanced Settings</string>
</property>
</widget>
</item>
<item>
<widget class="QFrame" name="frame_2">
<property name="frameShape">
<enum>QFrame::StyledPanel</enum>
</property>
<property name="frameShadow">
<enum>QFrame::Raised</enum>
</property>
<layout class="QFormLayout" name="formLayout">
<property name="fieldGrowthPolicy">
<enum>QFormLayout::AllNonFixedFieldsGrow</enum>
</property>
<item row="0" column="0">
<widget class="QLabel" name="m_QBallReconstructionThresholdLabel_2">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>B0 Threshold</string>
</property>
<property name="wordWrap">
<bool>false</bool>
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QLineEdit" name="m_QBallReconstructionThreasholdEdit">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>0</string>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QCheckBox" name="m_OutputB0Image">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>Output B0-Image</string>
</property>
</widget>
</item>
<item row="3" column="0">
<widget class="QLabel" name="label_2">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>Spherical Harmonics:</string>
</property>
</widget>
</item>
<item row="5" column="0">
<widget class="QLabel" name="m_QBallReconstructionMaxLLevelTextLabel_2">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>Maximum l-Level</string>
</property>
<property name="wordWrap">
<bool>false</bool>
</property>
</widget>
</item>
<item row="5" column="1">
<widget class="QComboBox" name="m_QBallReconstructionMaxLLevelComboBox">
<property name="enabled">
<bool>true</bool>
</property>
<property name="currentIndex">
<number>-1</number>
</property>
</widget>
</item>
<item row="6" column="0">
<widget class="QLabel" name="m_QBallReconstructionLambdaTextLabel_2">
<property name="enabled">
<bool>true</bool>
</property>
<property name="toolTip">
<string>Regularization Parameter</string>
</property>
<property name="text">
<string>Lambda</string>
</property>
<property name="wordWrap">
<bool>false</bool>
</property>
</widget>
</item>
<item row="6" column="1">
<widget class="QLineEdit" name="m_QBallReconstructionLambdaLineEdit">
<property name="enabled">
<bool>true</bool>
</property>
<property name="text">
<string>0.006</string>
</property>
</widget>
</item>
<item row="4" column="0">
<widget class="QCheckBox" name="m_OutputCoeffsImage">
<property name="text">
<string>Output SH-Coeffs-Image</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<widget class="QComboBox" name="m_QBallReconstructionMethodComboBox">
<property name="currentIndex">
- <number>0</number>
+ <number>2</number>
</property>
<item>
<property name="text">
<string>Numerical</string>
</property>
</item>
<item>
<property name="text">
<string>Standard (SH)</string>
</property>
</item>
<item>
<property name="text">
<string>Solid Angle (SH)</string>
</property>
</item>
<item>
<property name="text">
<string>Constraint Solid Angle (SH)</string>
</property>
</item>
<item>
<property name="text">
<string>ADC-Profile only (SH)</string>
</property>
</item>
<item>
<property name="text">
<string>Raw Signal only (SH)</string>
</property>
</item>
<item>
<property name="text">
- <string>Mulit q-Ball (SH)</string>
+ <string>Multi-Shell (SH)</string>
</property>
</item>
</widget>
</item>
<item>
<widget class="QLabel" name="m_Description">
<property name="text">
<string>TextLabel</string>
</property>
</widget>
</item>
<item>
<widget class="QCommandLinkButton" name="m_ButtonStandard">
<property name="enabled">
<bool>false</bool>
</property>
<property name="toolTip">
<string/>
</property>
<property name="statusTip">
<string/>
</property>
<property name="whatsThis">
<string notr="true"/>
</property>
<property name="text">
<string>Start Reconstruction</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item>
<widget class="QGroupBox" name="m_QBallSelectionBox">
<property name="enabled">
<bool>true</bool>
</property>
<property name="layoutDirection">
<enum>Qt::LeftToRight</enum>
</property>
<property name="autoFillBackground">
<bool>false</bool>
</property>
<property name="title">
- <string>Multi q-Ball reconstruction</string>
+ <string>Multi-Shell Reconstruction</string>
</property>
<layout class="QVBoxLayout" name="verticalLayout_2"/>
</widget>
</item>
<item>
<spacer name="verticalSpacer">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>0</height>
</size>
</property>
</spacer>
</item>
</layout>
</widget>
<layoutdefault spacing="6" margin="11"/>
<resources/>
<connections/>
</ui>
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 <berryISelectionService.h>
#include <berryIWorkbenchWindow.h>
#include <berryIStructuredSelection.h>
// Qmitk
#include "QmitkStreamlineTrackingView.h"
#include "QmitkStdMultiWidget.h"
// Qt
#include <QMessageBox>
// MITK
#include <mitkImageToItk.h>
#include <mitkFiberBundleX.h>
#include <mitkImageCast.h>
// VTK
#include <vtkPolyData.h>
#include <vtkPoints.h>
#include <vtkCellArray.h>
#include <vtkSmartPointer.h>
#include <vtkPolyLine.h>
#include <vtkCellData.h>
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<mitk::DataNode*> 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<mitk::DataNode*>::iterator it = nodes.begin(); it != nodes.end(); ++it )
{
mitk::DataNode::Pointer node = *it;
if( node.IsNotNull() && dynamic_cast<mitk::Image*>(node->GetData()) )
{
if( dynamic_cast<mitk::TensorImage*>(node->GetData()) )
{
m_TensorImageNode = node;
m_TensorImage = dynamic_cast<mitk::TensorImage*>(node->GetData());
m_Controls->m_TensorImageLabel->setText(node->GetName().c_str());
}
else
{
bool isBinary = false;
node->GetPropertyValue<bool>("binary", isBinary);
if (isBinary && m_SeedRoi.IsNull())
{
m_SeedRoi = dynamic_cast<mitk::Image*>(node->GetData());
m_Controls->m_RoiImageLabel->setText(node->GetName().c_str());
}
else if (isBinary)
{
m_MaskImage = dynamic_cast<mitk::Image*>(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<float>, 3> TensorImageType;
typedef mitk::ImageToItk<TensorImageType> CastType;
typedef mitk::ImageToItk<ItkUCharImageType> 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<ItkUCharImageType>(m_SeedRoi, mask);
filter->SetSeedImage(mask);
}
if (m_MaskImage.IsNotNull())
{
ItkUCharImageType::Pointer mask = ItkUCharImageType::New();
mitk::CastToItkImage<ItkUCharImageType>(m_MaskImage, mask);
filter->SetMaskImage(mask);
}
filter->Update();
vtkSmartPointer<vtkPolyData> 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 @@
<?xml version="1.0" encoding="UTF-8"?>
<ui version="4.0">
<class>QmitkStreamlineTrackingViewControls</class>
<widget class="QWidget" name="QmitkStreamlineTrackingViewControls">
<property name="geometry">
<rect>
<x>0</x>
<y>0</y>
<width>480</width>
<height>553</height>
</rect>
</property>
<property name="minimumSize">
<size>
<width>0</width>
<height>0</height>
</size>
</property>
<property name="windowTitle">
<string>QmitkTemplate</string>
</property>
<layout class="QGridLayout" name="gridLayout_3">
<property name="topMargin">
<number>3</number>
</property>
<property name="bottomMargin">
<number>3</number>
</property>
<property name="spacing">
<number>0</number>
</property>
<item row="6" column="0">
<spacer name="spacer1">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeType">
<enum>QSizePolicy::Expanding</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>220</height>
</size>
</property>
</spacer>
</item>
<item row="3" column="0">
<widget class="QGroupBox" name="groupBox_2">
<property name="title">
<string>Parameters</string>
</property>
<layout class="QGridLayout" name="gridLayout_2">
<item row="6" column="1">
<widget class="QSlider" name="m_SeedsPerVoxelSlider">
<property name="toolTip">
<string>Number of tracts started in each voxel of the seed ROI.</string>
</property>
<property name="minimum">
<number>1</number>
</property>
<property name="maximum">
- <number>10</number>
+ <number>100</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="QSlider" name="m_AngularThresholdSlider">
<property name="toolTip">
<string>Fractional Anisotropy Threshold</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
<number>90</number>
</property>
<property name="value">
<number>45</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="5" column="1">
<widget class="QSlider" name="m_MinTractLengthSlider">
<property name="toolTip">
<string>Minimum tract length in mm.</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
<number>500</number>
</property>
<property name="value">
<number>40</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="6" column="0">
<widget class="QLabel" name="m_SeedsPerVoxelLabel">
<property name="toolTip">
<string>Number of tracts started in each voxel of the seed ROI.</string>
</property>
<property name="text">
<string>Seeds per Voxel: 1</string>
</property>
</widget>
</item>
<item row="7" column="0">
<spacer name="horizontalSpacer">
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
<property name="sizeType">
<enum>QSizePolicy::Fixed</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>200</width>
<height>0</height>
</size>
</property>
</spacer>
</item>
<item row="0" column="1">
<widget class="QSlider" name="m_FaThresholdSlider">
<property name="toolTip">
<string>Fractional Anisotropy Threshold</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
<number>100</number>
</property>
<property name="value">
<number>20</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="2" column="1">
<widget class="QSlider" name="m_fSlider">
<property name="toolTip">
<string>Weighting factor between first eigenvector (f=1 equals FACT tracking) and input vector dependent direction (f=0).</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
<number>100</number>
</property>
<property name="value">
<number>100</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="0" column="0">
<widget class="QLabel" name="m_FaThresholdLabel">
<property name="toolTip">
<string>Minimum tract length in mm.</string>
</property>
<property name="text">
<string>FA Threshold: 0.2</string>
</property>
</widget>
</item>
<item row="5" column="0">
<widget class="QLabel" name="m_MinTractLengthLabel">
<property name="toolTip">
<string>Minimum tract length in mm.</string>
</property>
<property name="text">
<string>Min. Tract Length: 40mm</string>
</property>
</widget>
</item>
<item row="4" column="1">
<widget class="QSlider" name="m_StepsizeSlider">
<property name="toolTip">
<string>Stepsize in mm (auto = 0.1*minimal spacing)</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
<number>100</number>
</property>
<property name="value">
<number>0</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="3" column="0">
<widget class="QLabel" name="m_gLabel">
<property name="toolTip">
<string>Minimum tract length in mm.</string>
</property>
<property name="text">
<string>g: 0</string>
</property>
</widget>
</item>
<item row="3" column="1">
<widget class="QSlider" name="m_gSlider">
<property name="toolTip">
<string>Weighting factor between input vector (g=0) and tensor deflection (g=1 equals TEND tracking)</string>
</property>
<property name="minimum">
<number>0</number>
</property>
<property name="maximum">
<number>100</number>
</property>
<property name="value">
<number>0</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="m_AngularThresholdLabel">
<property name="toolTip">
<string>Minimum tract length in mm.</string>
</property>
<property name="text">
<string>Angular Threshold: 45°</string>
</property>
</widget>
</item>
<item row="4" column="0">
<widget class="QLabel" name="m_StepsizeLabel">
<property name="toolTip">
<string>Minimum tract length in mm.</string>
</property>
<property name="text">
<string>Step Size: auto</string>
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QLabel" name="m_fLabel">
<property name="toolTip">
<string>Minimum tract length in mm.</string>
</property>
<property name="text">
<string>f: 1</string>
</property>
</widget>
</item>
<item row="8" column="0">
<widget class="QCheckBox" name="m_InterpolationBox">
<property name="toolTip">
<string>Default is nearest neighbor interpolation.</string>
</property>
<property name="text">
<string>Enable trilinear interpolation</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item row="7" column="0">
<widget class="QLabel" name="m_SeedsPerVoxelLabel_2">
<property name="toolTip">
<string>Number of tracts started in each voxel of the seed ROI.</string>
</property>
<property name="text">
<string>Mori et al. Annals Neurology 1999</string>
</property>
</widget>
</item>
<item row="5" column="0">
<widget class="QCommandLinkButton" name="commandLinkButton">
<property name="enabled">
<bool>false</bool>
</property>
<property name="text">
<string>Start Tracking</string>
</property>
</widget>
</item>
<item row="0" column="0">
<widget class="QGroupBox" name="groupBox">
<property name="title">
<string>Data</string>
</property>
<layout class="QGridLayout" name="gridLayout">
<item row="0" column="1">
<widget class="QLabel" name="m_TensorImageLabel">
<property name="text">
<string>-</string>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="QLabel" name="m_RoiImageLabel">
<property name="text">
<string>-</string>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="label_6">
<property name="text">
<string>Seed ROI Image:</string>
</property>
</widget>
</item>
<item row="0" column="0">
<widget class="QLabel" name="label_2">
<property name="text">
<string>Tensor Image:</string>
</property>
</widget>
</item>
<item row="2" column="0">
<widget class="QLabel" name="label_7">
<property name="toolTip">
<string>Stop tracking if leaving mask</string>
</property>
<property name="text">
<string>Mask Image:</string>
</property>
</widget>
</item>
<item row="2" column="1">
<widget class="QLabel" name="m_MaskImageLabel">
<property name="text">
<string>-</string>
</property>
</widget>
</item>
</layout>
</widget>
</item>
<item row="9" column="0">
<widget class="QLabel" name="m_SeedsPerVoxelLabel_3">
<property name="toolTip">
<string>Number of tracts started in each voxel of the seed ROI.</string>
</property>
<property name="text">
<string>Lazar et al. Human Brain Mapping 2003</string>
</property>
</widget>
</item>
<item row="8" column="0">
<widget class="QLabel" name="m_SeedsPerVoxelLabel_4">
<property name="toolTip">
<string>Number of tracts started in each voxel of the seed ROI.</string>
</property>
<property name="text">
<string>Weinstein et al. Proceedings of IEEE Visualization 1999</string>
</property>
</widget>
</item>
</layout>
</widget>
<layoutdefault spacing="6" margin="11"/>
<tabstops>
<tabstop>commandLinkButton</tabstop>
</tabstops>
<resources/>
<connections/>
</ui>