diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp index d3c70e5165..954a726f06 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp @@ -1,465 +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(ItkQBallImgType* qballImage, ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen) +EnergyComputer::EnergyComputer(ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen) : m_UseTrilinearInterpolation(true) { m_ParticleGrid = particleGrid; m_RandGen = randGen; - m_Image = qballImage; + //m_Image = qballImage; m_SphereInterpolator = interpolator; m_Mask = mask; m_ParticleLength = m_ParticleGrid->m_ParticleLength; m_SquaredParticleLength = m_ParticleLength*m_ParticleLength; - m_Size[0] = m_Image->GetLargestPossibleRegion().GetSize()[0]; - m_Size[1] = m_Image->GetLargestPossibleRegion().GetSize()[1]; - m_Size[2] = m_Image->GetLargestPossibleRegion().GetSize()[2]; + 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] = m_Image->GetSpacing()[0]; - m_Spacing[1] = m_Image->GetSpacing()[1]; - m_Spacing[2] = m_Image->GetSpacing()[2]; + m_Spacing[0] = mask->GetSpacing()[0]; + m_Spacing[1] = mask->GetSpacing()[1]; + m_Spacing[2] = mask->GetSpacing()[2]; // calculate rotation matrix - vnl_matrix temp = m_Image->GetDirection().GetVnlMatrix(); + 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"); 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::EvaluateOdf(vnl_vector_fixed& pos, vnl_vector_fixed dir) -{ - const int sampleSteps = 10; // evaluate ODF at 2*sampleSteps+1 positions along dir - vnl_vector_fixed samplePos; // current position to evaluate - float result = 0; // average of sampled ODF values - int xint, yint, zint; // voxel containing samplePos - - // rotate particle direction according to image rotation - dir = m_RotationMatrix*dir; - - // get interpolation for rotated direction - m_SphereInterpolator->getInterpolation(dir); - - // sample ODF values along particle direction - for (int i=-sampleSteps; i <= sampleSteps;i++) - { - samplePos = pos + (dir * m_ParticleLength) * ((float)i/sampleSteps); - - if (!m_UseTrilinearInterpolation) // image has not enough slices to use trilinear interpolation - { - ItkQBallImgType::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_Image->GetLargestPossibleRegion().IsInside(index)) - { - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2]); - } - } - else // use trilinear interpolation - { - float Rx = samplePos[0]/m_Spacing[0]-0.5; - float Ry = samplePos[1]/m_Spacing[1]-0.5; - float Rz = samplePos[2]/m_Spacing[2]-0.5; - - xint = floor(Rx); - yint = floor(Ry); - zint = floor(Rz); - - if (xint >= 0 && xint < m_Size[0]-1 && yint >= 0 && yint < m_Size[1]-1 && zint >= 0 && zint < m_Size[2]-1) - { - float xfrac = Rx-xint; - float yfrac = Ry-yint; - float zfrac = Rz-zint; - - ItkQBallImgType::IndexType index; - float weight; - - weight = (1-xfrac)*(1-yfrac)*(1-zfrac); - index[0] = xint; index[1] = yint; index[2] = zint; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - - weight = (xfrac)*(1-yfrac)*(1-zfrac); - index[0] = xint+1; index[1] = yint; index[2] = zint; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - - weight = (1-xfrac)*(yfrac)*(1-zfrac); - index[0] = xint; index[1] = yint+1; index[2] = zint; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - - weight = (1-xfrac)*(1-yfrac)*(zfrac); - index[0] = xint; index[1] = yint; index[2] = zint+1; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - - weight = (xfrac)*(yfrac)*(1-zfrac); - index[0] = xint+1; index[1] = yint+1; index[2] = zint; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - - weight = (1-xfrac)*(yfrac)*(zfrac); - index[0] = xint; index[1] = yint+1; index[2] = zint+1; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - - weight = (xfrac)*(1-yfrac)*(zfrac); - index[0] = xint+1; index[1] = yint; index[2] = zint+1; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - - weight = (xfrac)*(yfrac)*(zfrac); - index[0] = xint+1; index[1] = yint+1; index[2] = zint+1; - result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; - } - } - } - result /= (2*sampleSteps+1); // average result over taken samples - return result; -} - -float EnergyComputer::ComputeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, Particle *dp) -{ - if (SpatProb(R) == 0) // check if position is inside mask - return itk::NumericTraits::NonpositiveMin(); - - float odfVal = EvaluateOdf(R, N); // evaluate ODF in given direction - - float modelVal = 0; - m_ParticleGrid->ComputeNeighbors(R); // retrieve neighbouring particles from particle grid - Particle* neighbour = m_ParticleGrid->GetNextNeighbor(); - while (neighbour!=NULL) // iterate over nieghbouring particles - { - if (dp != neighbour) // don't evaluate against itself - { - // see Reisert et al. "Global Reconstruction of Neuronal Fibers", MICCAI 2009 - float dot = fabs(dot_product(N,neighbour->dir)); - float bw = mbesseli0(dot); - float dpos = (neighbour->pos-R).squared_magnitude(); - float w = mexp(dpos*gamma_s); - modelVal += w*(bw+m_ParticleChemicalPotential); - w = mexp(dpos*gamma_reg_s); - } - neighbour = m_ParticleGrid->GetNextNeighbor(); - } - - float energy = 2*(odfVal/m_ParticleWeight-modelVal) - (mbesseli0(1.0)+m_ParticleChemicalPotential); - return energy*m_ExtStrength; -} - -float EnergyComputer::ComputeInternalEnergy(Particle *dp) -{ - float energy = 0; - - if (dp->pID != -1) // has predecessor - energy += ComputeInternalEnergyConnection(dp,+1); - - if (dp->mID != -1) // has successor - energy += ComputeInternalEnergyConnection(dp,-1); - - return energy; -} - -float EnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1) -{ - Particle *p2 = 0; - int ep2; - - if (ep1 == 1) - p2 = m_ParticleGrid->GetParticle(p1->pID); // get predecessor - else - p2 = m_ParticleGrid->GetParticle(p1->mID); // get successor - - // check in which direction the connected particle is pointing - if (p2->mID == p1->ID) - ep2 = -1; - else if (p2->pID == p1->ID) - ep2 = 1; - else - std::cout << "EnergyComputer: Connections are inconsistent!" << std::endl; - - return ComputeInternalEnergyConnection(p1,ep1,p2,ep2); -} - -float EnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) -{ - // see Reisert et al. "Global Reconstruction of Neuronal Fibers", MICCAI 2009 - if ((dot_product(p1->dir,p2->dir))*ep1*ep2 > -m_CurvatureThreshold) // angle between particles is too sharp - return itk::NumericTraits::NonpositiveMin(); - - // calculate the endpoints of the two particles - vnl_vector_fixed endPoint1 = p1->pos + (p1->dir * (m_ParticleLength * ep1)); - vnl_vector_fixed endPoint2 = p2->pos + (p2->dir * (m_ParticleLength * ep2)); - - // check if endpoints are too far apart to connect - if ((endPoint1-endPoint2).squared_magnitude() > m_SquaredParticleLength) - return itk::NumericTraits::NonpositiveMin(); - - // calculate center point of the two particles - vnl_vector_fixed R = (p2->pos + p1->pos); R *= 0.5; - - // they are not allowed to connect if the mask image does not allow it - if (SpatProb(R) == 0) - return itk::NumericTraits::NonpositiveMin(); - - // get distances of endpoints to center point - float norm1 = (endPoint1-R).squared_magnitude(); - float norm2 = (endPoint2-R).squared_magnitude(); - - // calculate actual internal energy - float energy = (m_ConnectionPotential-norm1-norm2)*m_IntStrength; - return energy; -} - 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; } - -//ComputeFiberCorrelation() -//{ -// float bD = 15; - -// vnl_matrix_fixed bDir = -// *itk::PointShell >::DistributePointShell(); - -// const int N = QBALL_ODFSIZE; - -// vnl_matrix_fixed temp = bDir.transpose(); -// vnl_matrix_fixed C = temp*bDir; -// vnl_matrix_fixed Q = C; -// vnl_vector_fixed mean; -// for(int i=0; i repMean; -// for (int i=0; i P = Q*Q; - -// std::vector pointer; -// pointer.reserve(N*N); -// double * start = C.data_block(); -// double * end = start + N*N; -// for (double * iter = start; iter != end; ++iter) -// { -// pointer.push_back(iter); -// } -// std::sort(pointer.begin(), pointer.end(), LessDereference()); - -// vnl_vector_fixed alpha; -// vnl_vector_fixed beta; -// for (int i=0; im_Meanval_sq = (sum*sum)/N; - -// vnl_vector_fixed alpha_0; -// vnl_vector_fixed alpha_2; -// vnl_vector_fixed alpha_4; -// vnl_vector_fixed alpha_6; -// for(int i=0; i T; -// T.set_column(0,alpha_0); -// T.set_column(1,alpha_2); -// T.set_column(2,alpha_4); -// T.set_column(3,alpha_6); - -// vnl_vector_fixed coeff = vnl_matrix_inverse(T).pinverse()*beta; - -// MITK_INFO << "itkGibbsTrackingFilter: Bessel oefficients: " << coeff; - -// BESSEL_APPROXCOEFF = new float[4]; - -// BESSEL_APPROXCOEFF[0] = coeff(0); -// BESSEL_APPROXCOEFF[1] = coeff(1); -// BESSEL_APPROXCOEFF[2] = coeff(2); -// BESSEL_APPROXCOEFF[3] = coeff(3); -// BESSEL_APPROXCOEFF[0] = -0.1714; -// BESSEL_APPROXCOEFF[1] = 0.5332; -// BESSEL_APPROXCOEFF[2] = -1.4889; -// BESSEL_APPROXCOEFF[3] = 2.0389; -//} diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h index f8b1f357af..d7452eff70 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h @@ -1,86 +1,86 @@ /*=================================================================== 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 _ENCOMP #define _ENCOMP #include #include #include #include #include using namespace mitk; class MitkDiffusionImaging_EXPORT EnergyComputer { public: - typedef itk::Vector OdfVectorType; - typedef itk::Image ItkQBallImgType; + //typedef itk::Vector OdfVectorType; + //typedef itk::Image ItkQBallImgType; typedef itk::Image ItkFloatImageType; typedef itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRandGenType; - EnergyComputer(ItkQBallImgType* qballImage, ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen); + EnergyComputer(ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen); void SetParameters(float particleWeight, float particleWidth, float connectionPotential, float curvThres, float inexBalance, float particlePotential); // get random position inside mask void DrawRandomPosition(vnl_vector_fixed& R); // external energy calculation - float ComputeExternalEnergy(vnl_vector_fixed& R, vnl_vector_fixed& N, Particle* dp); + virtual float ComputeExternalEnergy(vnl_vector_fixed& R, vnl_vector_fixed& N, Particle* dp) =0; // internal energy calculation - float ComputeInternalEnergyConnection(Particle *p1,int ep1); - float ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2); - float ComputeInternalEnergy(Particle *dp); + virtual float ComputeInternalEnergyConnection(Particle *p1,int ep1) = 0; + virtual float ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) = 0; + virtual float ComputeInternalEnergy(Particle *dp) = 0; int GetNumActiveVoxels(); protected: vnl_matrix_fixed m_RotationMatrix; SphereInterpolator* m_SphereInterpolator; ParticleGrid* m_ParticleGrid; ItkRandGenType* m_RandGen; - ItkQBallImgType* m_Image; +// ItkQBallImgType* m_Image; ItkFloatImageType* m_Mask; vnl_vector_fixed m_Size; vnl_vector_fixed m_Spacing; std::vector< float > m_CumulatedSpatialProbability; std::vector< int > m_ActiveIndices; // indices inside mask bool m_UseTrilinearInterpolation; // is deactivated if less than 3 image slices are available int m_NumActiveVoxels; // voxels inside mask float m_ConnectionPotential; // larger value results in larger energy value -> higher proposal acceptance probability float m_ParticleChemicalPotential; // larger value results in larger energy value -> higher proposal acceptance probability float gamma_s; float gamma_reg_s; float m_ParticleWeight; // defines how much one particle contributes to the artificial signal float m_ExtStrength; // weighting factor for external energy float m_IntStrength; // weighting factor for internal energy float m_ParticleLength; // particle length float m_SquaredParticleLength; // squared particle length float m_CurvatureThreshold; // maximum angle accepted between two connected particles float SpatProb(vnl_vector_fixed pos); float EvaluateOdf(vnl_vector_fixed &pos, vnl_vector_fixed dir); float mbesseli0(float x); float mexp(float x); }; #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkGibbsEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkGibbsEnergyComputer.cpp new file mode 100644 index 0000000000..8310e5cf48 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkGibbsEnergyComputer.cpp @@ -0,0 +1,213 @@ +#include + +#include +#include +#include +#include +#include + +using namespace mitk; + +GibbsEnergyComputer::GibbsEnergyComputer(ItkQBallImgType* qballImage, ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen) +:EnergyComputer(mask, particleGrid, interpolator, randGen) + +{ + m_Image = qballImage; +} + +float GibbsEnergyComputer::EvaluateOdf(vnl_vector_fixed& pos, vnl_vector_fixed dir) +{ + const int sampleSteps = 10; // evaluate ODF at 2*sampleSteps+1 positions along dir + vnl_vector_fixed samplePos; // current position to evaluate + float result = 0; // average of sampled ODF values + int xint, yint, zint; // voxel containing samplePos + + // rotate particle direction according to image rotation + dir = m_RotationMatrix*dir; + + // get interpolation for rotated direction + m_SphereInterpolator->getInterpolation(dir); + + // sample ODF values along particle direction + for (int i=-sampleSteps; i <= sampleSteps;i++) + { + samplePos = pos + (dir * m_ParticleLength) * ((float)i/sampleSteps); + + if (!m_UseTrilinearInterpolation) // image has not enough slices to use trilinear interpolation + { + ItkQBallImgType::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_Image->GetLargestPossibleRegion().IsInside(index)) + { + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2]); + } + } + else // use trilinear interpolation + { + float Rx = samplePos[0]/m_Spacing[0]-0.5; + float Ry = samplePos[1]/m_Spacing[1]-0.5; + float Rz = samplePos[2]/m_Spacing[2]-0.5; + + xint = floor(Rx); + yint = floor(Ry); + zint = floor(Rz); + + if (xint >= 0 && xint < m_Size[0]-1 && yint >= 0 && yint < m_Size[1]-1 && zint >= 0 && zint < m_Size[2]-1) + { + float xfrac = Rx-xint; + float yfrac = Ry-yint; + float zfrac = Rz-zint; + + ItkQBallImgType::IndexType index; + float weight; + + weight = (1-xfrac)*(1-yfrac)*(1-zfrac); + index[0] = xint; index[1] = yint; index[2] = zint; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + + weight = (xfrac)*(1-yfrac)*(1-zfrac); + index[0] = xint+1; index[1] = yint; index[2] = zint; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + + weight = (1-xfrac)*(yfrac)*(1-zfrac); + index[0] = xint; index[1] = yint+1; index[2] = zint; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + + weight = (1-xfrac)*(1-yfrac)*(zfrac); + index[0] = xint; index[1] = yint; index[2] = zint+1; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + + weight = (xfrac)*(yfrac)*(1-zfrac); + index[0] = xint+1; index[1] = yint+1; index[2] = zint; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + + weight = (1-xfrac)*(yfrac)*(zfrac); + index[0] = xint; index[1] = yint+1; index[2] = zint+1; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + + weight = (xfrac)*(1-yfrac)*(zfrac); + index[0] = xint+1; index[1] = yint; index[2] = zint+1; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + + weight = (xfrac)*(yfrac)*(zfrac); + index[0] = xint+1; index[1] = yint+1; index[2] = zint+1; + result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + } + } + } + result /= (2*sampleSteps+1); // average result over taken samples + return result; +} + +float GibbsEnergyComputer::ComputeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, Particle *dp) +{ + if (SpatProb(R) == 0) // check if position is inside mask + return itk::NumericTraits::NonpositiveMin(); + + float odfVal = EvaluateOdf(R, N); // evaluate ODF in given direction + + float modelVal = 0; + m_ParticleGrid->ComputeNeighbors(R); // retrieve neighbouring particles from particle grid + Particle* neighbour = m_ParticleGrid->GetNextNeighbor(); + while (neighbour!=NULL) // iterate over nieghbouring particles + { + if (dp != neighbour) // don't evaluate against itself + { + // see Reisert et al. "Global Reconstruction of Neuronal Fibers", MICCAI 2009 + float dot = fabs(dot_product(N,neighbour->dir)); + float bw = mbesseli0(dot); + float dpos = (neighbour->pos-R).squared_magnitude(); + float w = mexp(dpos*gamma_s); + modelVal += w*(bw+m_ParticleChemicalPotential); + w = mexp(dpos*gamma_reg_s); + } + neighbour = m_ParticleGrid->GetNextNeighbor(); + } + + float energy = 2*(odfVal/m_ParticleWeight-modelVal) - (mbesseli0(1.0)+m_ParticleChemicalPotential); + return energy*m_ExtStrength; +} + +float GibbsEnergyComputer::ComputeInternalEnergy(Particle *dp) +{ + float energy = 0; + + if (dp->pID != -1) // has predecessor + energy += ComputeInternalEnergyConnection(dp,+1); + + if (dp->mID != -1) // has successor + energy += ComputeInternalEnergyConnection(dp,-1); + + return energy; +} + +float GibbsEnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1) +{ + Particle *p2 = 0; + int ep2; + + if (ep1 == 1) + p2 = m_ParticleGrid->GetParticle(p1->pID); // get predecessor + else + p2 = m_ParticleGrid->GetParticle(p1->mID); // get successor + + // check in which direction the connected particle is pointing + if (p2->mID == p1->ID) + ep2 = -1; + else if (p2->pID == p1->ID) + ep2 = 1; + else + std::cout << "EnergyComputer: Connections are inconsistent!" << std::endl; + + return ComputeInternalEnergyConnection(p1,ep1,p2,ep2); +} + +float GibbsEnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) +{ + // see Reisert et al. "Global Reconstruction of Neuronal Fibers", MICCAI 2009 + if ((dot_product(p1->dir,p2->dir))*ep1*ep2 > -m_CurvatureThreshold) // angle between particles is too sharp + return itk::NumericTraits::NonpositiveMin(); + + // calculate the endpoints of the two particles + vnl_vector_fixed endPoint1 = p1->pos + (p1->dir * (m_ParticleLength * ep1)); + vnl_vector_fixed endPoint2 = p2->pos + (p2->dir * (m_ParticleLength * ep2)); + + // check if endpoints are too far apart to connect + if ((endPoint1-endPoint2).squared_magnitude() > m_SquaredParticleLength) + return itk::NumericTraits::NonpositiveMin(); + + // calculate center point of the two particles + vnl_vector_fixed R = (p2->pos + p1->pos); R *= 0.5; + + // they are not allowed to connect if the mask image does not allow it + if (SpatProb(R) == 0) + return itk::NumericTraits::NonpositiveMin(); + + // get distances of endpoints to center point + float norm1 = (endPoint1-R).squared_magnitude(); + float norm2 = (endPoint2-R).squared_magnitude(); + + // calculate actual internal energy + float energy = (m_ConnectionPotential-norm1-norm2)*m_IntStrength; + return energy; +} \ No newline at end of file diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkGibbsEnergyComputer.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkGibbsEnergyComputer.h new file mode 100644 index 0000000000..c2d462042d --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkGibbsEnergyComputer.h @@ -0,0 +1,41 @@ +#ifndef GIBBSENERGYCOMPUTER_H +#define GIBBSENERGYCOMPUTER_H + +#include +#include +#include +#include +#include + +#include "mitkEnergyComputer.h" + + +using namespace mitk; + +class MitkDiffusionImaging_EXPORT GibbsEnergyComputer : public EnergyComputer +{ + public: + + typedef itk::Vector OdfVectorType; + typedef itk::Image ItkQBallImgType; + typedef itk::Image ItkFloatImageType; + typedef itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRandGenType; + + GibbsEnergyComputer(ItkQBallImgType* qballImage, ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen); + + // external energy calculation + float ComputeExternalEnergy(vnl_vector_fixed& R, vnl_vector_fixed& N, Particle* dp); + + // internal energy calculation + float ComputeInternalEnergyConnection(Particle *p1,int ep1); + float ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2); + float ComputeInternalEnergy(Particle *dp); + + float EvaluateOdf(vnl_vector_fixed& pos, vnl_vector_fixed dir); + protected: + + ItkQBallImgType* m_Image; + +}; + +#endif \ No newline at end of file diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index 3c93e73ca5..95c38da0e0 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,411 +1,416 @@ /*=================================================================== 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 "itkGibbsTrackingFilter.h" // MITK #include #include #include #include #include -#include +//#include #include +#include // ITK #include #include // MISC #include #include #include #include namespace itk{ template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::GibbsTrackingFilter(): m_StartTemperature(0.1), m_EndTemperature(0.001), m_Iterations(500000), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), m_ConnectionPotential(10), m_InexBalance(0), m_ParticlePotential(0.2), m_MinFiberLength(10), m_AbortTracking(false), m_NumConnections(0), m_NumParticles(0), m_NumAcceptedFibers(0), m_CurrentStep(0), m_BuildFibers(false), m_Steps(10), m_ProposalAcceptance(0), m_CurvatureThreshold(0.7), m_DuplicateImage(true), m_RandomSeed(-1), m_ParameterFile(""), m_LutPath("") { } template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::~GibbsTrackingFilter() { } // fill output fiber bundle datastructure template< class ItkQBallImageType > typename GibbsTrackingFilter< ItkQBallImageType >::FiberPolyDataType GibbsTrackingFilter< ItkQBallImageType >::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } return m_FiberPolyData; } template< class ItkQBallImageType > void GibbsTrackingFilter< ItkQBallImageType > ::EstimateParticleWeight() { MITK_INFO << "GibbsTrackingFilter: estimating particle weight"; float minSpacing; if(m_QBallImage->GetSpacing()[0]GetSpacing()[1] && m_QBallImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[0]; else if (m_QBallImage->GetSpacing()[1] < m_QBallImage->GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[1]; else minSpacing = m_QBallImage->GetSpacing()[2]; float m_ParticleLength = 1.5*minSpacing; float m_ParticleWidth = 0.5*minSpacing; // seed random generators Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = Statistics::MersenneTwisterRandomVariateGenerator::New(); if (m_RandomSeed>-1) randGen->SetSeed(m_RandomSeed); else randGen->SetSeed(); // instantiate all necessary components SphereInterpolator* interpolator = new SphereInterpolator(m_LutPath); ParticleGrid* particleGrid = new ParticleGrid(m_MaskImage, m_ParticleLength); - EnergyComputer* encomp = new EnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); + GibbsEnergyComputer* encomp = new GibbsEnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); + + // EnergyComputer* encomp = new EnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); MetropolisHastingsSampler* sampler = new MetropolisHastingsSampler(particleGrid, encomp, randGen, m_CurvatureThreshold); float alpha = log(m_EndTemperature/m_StartTemperature); m_ParticleWeight = 0.01; int ppv = 0; // main loop int neededParts = 3000; while (ppvSetParameters(m_ParticleWeight,m_ParticleWidth,m_ConnectionPotential*m_ParticleLength*m_ParticleLength,m_CurvatureThreshold,m_InexBalance,m_ParticlePotential); for( int step = 0; step < 10; step++ ) { // update temperatur for simulated annealing process float temperature = m_StartTemperature * exp(alpha*(((1.0)*step)/((1.0)*10))); sampler->SetTemperature(temperature); for (unsigned long i=0; i<10000; i++) sampler->MakeProposal(); } ppv = particleGrid->m_NumParticles; particleGrid->ResetGrid(); } delete sampler; delete encomp; delete particleGrid; delete interpolator; MITK_INFO << "GibbsTrackingFilter: finished estimating particle weight"; } // perform global tracking template< class ItkQBallImageType > void GibbsTrackingFilter< ItkQBallImageType >::GenerateData() { // check if input is qball or tensor image and generate qball if necessary if (m_QBallImage.IsNull() && m_TensorImage.IsNotNull()) { TensorImageToQBallImageFilter::Pointer filter = TensorImageToQBallImageFilter::New(); filter->SetInput( m_TensorImage ); filter->Update(); m_QBallImage = filter->GetOutput(); } else if (m_DuplicateImage) // generate local working copy of QBall image (if not disabled) { typedef itk::ImageDuplicator< ItkQBallImageType > DuplicateFilterType; typename DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); duplicator->SetInputImage( m_QBallImage ); duplicator->Update(); m_QBallImage = duplicator->GetOutput(); } // perform mean subtraction on odfs typedef ImageRegionIterator< ItkQBallImageType > InputIteratorType; InputIteratorType it(m_QBallImage, m_QBallImage->GetLargestPossibleRegion() ); it.GoToBegin(); while (!it.IsAtEnd()) { itk::OrientationDistributionFunction odf(it.Get().GetDataPointer()); float mean = odf.GetMeanValue(); odf -= mean; it.Set(odf.GetDataPointer()); ++it; } // check if mask image is given if it needs resampling PrepareMaskImage(); // load parameter file LoadParameters(m_ParameterFile); // prepare parameters float minSpacing; if(m_QBallImage->GetSpacing()[0]GetSpacing()[1] && m_QBallImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[0]; else if (m_QBallImage->GetSpacing()[1] < m_QBallImage->GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[1]; else minSpacing = m_QBallImage->GetSpacing()[2]; if(m_ParticleLength == 0) m_ParticleLength = 1.5*minSpacing; if(m_ParticleWidth == 0) m_ParticleWidth = 0.5*minSpacing; if(m_ParticleWeight == 0) EstimateParticleWeight(); float alpha = log(m_EndTemperature/m_StartTemperature); m_Steps = m_Iterations/10000; if (m_Steps<10) m_Steps = 10; if (m_Steps>m_Iterations) { MITK_INFO << "GibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } if (m_CurvatureThreshold < mitk::eps) m_CurvatureThreshold = 0; unsigned long singleIts = (unsigned long)((1.0*m_Iterations) / (1.0*m_Steps)); // seed random generators Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = Statistics::MersenneTwisterRandomVariateGenerator::New(); if (m_RandomSeed>-1) randGen->SetSeed(m_RandomSeed); else randGen->SetSeed(); // load sphere interpolator to evaluate the ODFs SphereInterpolator* interpolator = new SphereInterpolator(m_LutPath); // initialize the actual tracking components (ParticleGrid, Metropolis Hastings Sampler and Energy Computer) ParticleGrid* particleGrid = new ParticleGrid(m_MaskImage, m_ParticleLength); - EnergyComputer* encomp = new EnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); + GibbsEnergyComputer* encomp = new GibbsEnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); + +// EnergyComputer* encomp = new EnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); encomp->SetParameters(m_ParticleWeight,m_ParticleWidth,m_ConnectionPotential*m_ParticleLength*m_ParticleLength,m_CurvatureThreshold,m_InexBalance,m_ParticlePotential); MetropolisHastingsSampler* sampler = new MetropolisHastingsSampler(particleGrid, encomp, randGen, m_CurvatureThreshold); MITK_INFO << "----------------------------------------"; MITK_INFO << "Iterations: " << m_Iterations; MITK_INFO << "Steps: " << m_Steps; MITK_INFO << "Particle length: " << m_ParticleLength; MITK_INFO << "Particle width: " << m_ParticleWidth; MITK_INFO << "Particle weight: " << m_ParticleWeight; MITK_INFO << "Start temperature: " << m_StartTemperature; MITK_INFO << "End temperature: " << m_EndTemperature; MITK_INFO << "In/Ex balance: " << m_InexBalance; MITK_INFO << "Min. fiber length: " << m_MinFiberLength; MITK_INFO << "Curvature threshold: " << m_CurvatureThreshold; MITK_INFO << "Random seed: " << m_RandomSeed; MITK_INFO << "----------------------------------------"; // main loop m_NumAcceptedFibers = 0; unsigned long counter = 1; for( m_CurrentStep = 1; m_CurrentStep <= m_Steps; m_CurrentStep++ ) { // update temperatur for simulated annealing process float temperature = m_StartTemperature * exp(alpha*(((1.0)*m_CurrentStep)/((1.0)*m_Steps))); sampler->SetTemperature(temperature); for (unsigned long i=0; iMakeProposal(); if (m_BuildFibers || (i==singleIts-1 && m_CurrentStep==m_Steps)) { m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/counter; m_NumParticles = particleGrid->m_NumParticles; m_NumConnections = particleGrid->m_NumConnections; FiberBuilder fiberBuilder(particleGrid, m_MaskImage); m_FiberPolyData = fiberBuilder.iterate(m_MinFiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); m_BuildFibers = false; } counter++; } m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/counter; m_NumParticles = particleGrid->m_NumParticles; m_NumConnections = particleGrid->m_NumConnections; MITK_INFO << "GibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; MITK_INFO << "GibbsTrackingFilter: particles: " << m_NumParticles; MITK_INFO << "GibbsTrackingFilter: connections: " << m_NumConnections; MITK_INFO << "GibbsTrackingFilter: progress: " << 100*(float)m_CurrentStep/m_Steps << "%"; MITK_INFO << "GibbsTrackingFilter: cell overflows: " << particleGrid->m_NumCellOverflows; MITK_INFO << "----------------------------------------"; if (m_AbortTracking) break; } delete sampler; delete encomp; delete interpolator; delete particleGrid; m_AbortTracking = true; m_BuildFibers = false; MITK_INFO << "GibbsTrackingFilter: done generate data"; } template< class ItkQBallImageType > void GibbsTrackingFilter< ItkQBallImageType >::PrepareMaskImage() { if(m_MaskImage.IsNull()) { MITK_INFO << "GibbsTrackingFilter: generating default mask image"; m_MaskImage = ItkFloatImageType::New(); m_MaskImage->SetSpacing( m_QBallImage->GetSpacing() ); m_MaskImage->SetOrigin( m_QBallImage->GetOrigin() ); m_MaskImage->SetDirection( m_QBallImage->GetDirection() ); m_MaskImage->SetRegions( m_QBallImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1.0); } else if ( m_MaskImage->GetLargestPossibleRegion().GetSize()[0]!=m_QBallImage->GetLargestPossibleRegion().GetSize()[0] || m_MaskImage->GetLargestPossibleRegion().GetSize()[1]!=m_QBallImage->GetLargestPossibleRegion().GetSize()[1] || m_MaskImage->GetLargestPossibleRegion().GetSize()[2]!=m_QBallImage->GetLargestPossibleRegion().GetSize()[2] || m_MaskImage->GetSpacing()[0]!=m_QBallImage->GetSpacing()[0] || m_MaskImage->GetSpacing()[1]!=m_QBallImage->GetSpacing()[1] || m_MaskImage->GetSpacing()[2]!=m_QBallImage->GetSpacing()[2] ) { MITK_INFO << "GibbsTrackingFilter: resampling mask image"; typedef itk::ResampleImageFilter< ItkFloatImageType, ItkFloatImageType, float > ResamplerType; ResamplerType::Pointer resampler = ResamplerType::New(); resampler->SetOutputSpacing( m_QBallImage->GetSpacing() ); resampler->SetOutputOrigin( m_QBallImage->GetOrigin() ); resampler->SetOutputDirection( m_QBallImage->GetDirection() ); resampler->SetSize( m_QBallImage->GetLargestPossibleRegion().GetSize() ); resampler->SetInput( m_MaskImage ); resampler->SetDefaultPixelValue(1.0); resampler->Update(); m_MaskImage = resampler->GetOutput(); MITK_INFO << "GibbsTrackingFilter: resampling finished"; } } // load current tracking paramters from xml file (.gtp) template< class ItkQBallImageType > bool GibbsTrackingFilter< ItkQBallImageType >::LoadParameters(std::string filename) { m_AbortTracking = true; try { if( filename.length()==0 ) { m_AbortTracking = false; return true; } MITK_INFO << "GibbsTrackingFilter: loading parameter file " << filename; TiXmlDocument doc( filename ); doc.LoadFile(); TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("parameter_set").Element(); QString iterations(pElem->Attribute("iterations")); m_Iterations = iterations.toULong(); QString particleLength(pElem->Attribute("particle_length")); m_ParticleLength = particleLength.toFloat(); QString particleWidth(pElem->Attribute("particle_width")); m_ParticleWidth = particleWidth.toFloat(); QString partWeight(pElem->Attribute("particle_weight")); m_ParticleWeight = partWeight.toFloat(); QString startTemp(pElem->Attribute("temp_start")); m_StartTemperature = startTemp.toFloat(); QString endTemp(pElem->Attribute("temp_end")); m_EndTemperature = endTemp.toFloat(); QString inExBalance(pElem->Attribute("inexbalance")); m_InexBalance = inExBalance.toFloat(); QString fiberLength(pElem->Attribute("fiber_length")); m_MinFiberLength = fiberLength.toFloat(); QString curvThres(pElem->Attribute("curvature_threshold")); m_CurvatureThreshold = cos(curvThres.toFloat()*M_PI/180); m_AbortTracking = false; MITK_INFO << "GibbsTrackingFilter: parameter file loaded successfully"; return true; } catch(...) { MITK_INFO << "GibbsTrackingFilter: could not load parameter file"; return false; } } } diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 5912a39400..e1ba28e8cb 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,239 +1,241 @@ set(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures -> Tbss IODataStructures/TbssImages/mitkTbssImageSource.cpp IODataStructures/TbssImages/mitkTbssRoiImageSource.cpp IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageIOFactory.cpp IODataStructures/TbssImages/mitkTbssImage.cpp IODataStructures/TbssImages/mitkTbssRoiImage.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriterFactory.cpp IODataStructures/TbssImages/mitkTbssImporter.cpp # DataStructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp Rendering/mitkTbssImageMapper.cpp Rendering/mitkPlanarCircleMapper3D.cpp Rendering/mitkPlanarPolygonMapper3D.cpp Rendering/mitkConnectomicsNetworkMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/mitkTractAnalyzer.cpp # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.cpp Algorithms/Connectomics/mitkConnectomicsHistogramBase.cpp Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.cpp Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.cpp Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.cpp Algorithms/Connectomics/mitkConnectomicsHistogramCache.cpp Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.cpp # Tractography Tractography/GibbsTracking/mitkParticleGrid.cpp Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp Tractography/GibbsTracking/mitkEnergyComputer.cpp + Tractography/GibbsTracking/mitkGibbsEnergyComputer.cpp Tractography/GibbsTracking/mitkFiberBuilder.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkTbssImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h Rendering/mitkPlanarCircleMapper3D.h Rendering/mitkPlanarPolygonMapper3D.h Rendering/mitkConnectomicsNetworkMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h Reconstruction/itkRegularizedIVIMReconstructionFilter.h Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/TbssImages/mitkTbssImporter.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h # Datastructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.h # Tractography Tractography/itkGibbsTrackingFilter.h Tractography/itkStochasticTractographyFilter.h Tractography/itkStreamlineTrackingFilter.h Tractography/GibbsTracking/mitkParticle.h Tractography/GibbsTracking/mitkParticleGrid.h Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h Tractography/GibbsTracking/mitkSimpSamp.h Tractography/GibbsTracking/mitkEnergyComputer.h + Tractography/GibbsTracking/mitkGibbsEnergyComputer.h Tractography/GibbsTracking/mitkSphereInterpolator.h Tractography/GibbsTracking/mitkFiberBuilder.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkSkeletonizationFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h Algorithms/itkElectrostaticRepulsionDiffusionGradientReductionFilter.h # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.h Algorithms/Connectomics/mitkConnectomicsHistogramBase.h Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.h Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.h Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.h Algorithms/Connectomics/mitkConnectomicsHistogramCache.h Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.h ) set( TOOL_FILES ) if(WIN32) endif(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp )