diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.cpp index c4c8fb3764..35df1b203a 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.cpp @@ -1,172 +1,171 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkEnergyComputer.h" #include #include using namespace mitk; EnergyComputer::EnergyComputer(ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen) : m_UseTrilinearInterpolation(true) { m_ParticleGrid = particleGrid; m_RandGen = randGen; //m_Image = qballImage; m_SphereInterpolator = interpolator; m_Mask = mask; m_ParticleLength = m_ParticleGrid->m_ParticleLength; m_SquaredParticleLength = m_ParticleLength*m_ParticleLength; m_Size[0] = mask->GetLargestPossibleRegion().GetSize()[0]; m_Size[1] = mask->GetLargestPossibleRegion().GetSize()[1]; m_Size[2] = mask->GetLargestPossibleRegion().GetSize()[2]; if (m_Size[0]<3 || m_Size[1]<3 || m_Size[2]<3) m_UseTrilinearInterpolation = false; m_Spacing[0] = mask->GetSpacing()[0]; m_Spacing[1] = mask->GetSpacing()[1]; m_Spacing[2] = mask->GetSpacing()[2]; // get rotation matrix vnl_matrix temp = mask->GetDirection().GetVnlMatrix(); vnl_matrix directionMatrix; directionMatrix.set_size(3,3); vnl_copy(temp, directionMatrix); vnl_matrix_fixed I = directionMatrix*directionMatrix.transpose(); for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) I[i][j] = fabs(I[i][j]); if(!I.is_identity(0.001)) { fprintf(stderr,"EnergyComputer: image direction is not a rotation matrix. Tracking not possible!\n"); std::cout << directionMatrix << std::endl; } 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); + m_SigmaInv = 1/(2*sigma_s*sigma_s); } // draw random position from active voxels void EnergyComputer::DrawRandomPosition(vnl_vector_fixed& R) { float r = m_RandGen->GetVariate();//m_RandGen->frand(); int j; int rl = 1; int rh = m_NumActiveVoxels; while(rh != rl) { j = rl + (rh-rl)/2; if (r < m_CumulatedSpatialProbability[j]) { rh = j; continue; } if (r > m_CumulatedSpatialProbability[j]) { rl = j+1; continue; } break; } R[0] = m_Spacing[0]*((float)(m_ActiveIndices[rh-1] % m_Size[0]) + m_RandGen->GetVariate()); R[1] = m_Spacing[1]*((float)((m_ActiveIndices[rh-1]/m_Size[0]) % m_Size[1]) + m_RandGen->GetVariate()); R[2] = m_Spacing[2]*((float)(m_ActiveIndices[rh-1]/(m_Size[0]*m_Size[1])) + m_RandGen->GetVariate()); } // return spatial probability of position float EnergyComputer::SpatProb(vnl_vector_fixed pos) { ItkFloatImageType::IndexType index; index[0] = floor(pos[0]/m_Spacing[0]); index[1] = floor(pos[1]/m_Spacing[1]); index[2] = floor(pos[2]/m_Spacing[2]); if (m_Mask->GetLargestPossibleRegion().IsInside(index)) // is inside image? return m_Mask->GetPixel(index); else return 0; } -float EnergyComputer::mbesseli0(float x) +float EnergyComputer::CalcI0(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/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.h index 88ade0c0bf..50d92b24a8 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkEnergyComputer.h @@ -1,88 +1,87 @@ /*=================================================================== 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; /** * \brief Calculates internal and external energy of the new particle configuration proposal. */ class MITKFIBERTRACKING_EXPORT EnergyComputer { public: //typedef itk::Vector OdfVectorType; //typedef itk::Image ItkQBallImgType; typedef itk::Image ItkFloatImageType; typedef itk::Statistics::MersenneTwisterRandomVariateGenerator ItkRandGenType; 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 virtual float ComputeExternalEnergy(vnl_vector_fixed& R, vnl_vector_fixed& N, Particle* dp) =0; // internal energy calculation 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; 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_SigmaInv; 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 CalcI0(float x); float mexp(float x); }; #endif diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp index 4cd2fcc8fa..48cf80a0f8 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp @@ -1,217 +1,219 @@ #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; } GibbsEnergyComputer::~GibbsEnergyComputer() { } 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!=nullptr) // 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->GetDir())); - float bw = mbesseli0(dot); + float dot = fabs( dot_product(N,neighbour->GetDir()) ); + float I_0 = CalcI0(dot) + m_ParticleChemicalPotential; + float dpos = (neighbour->GetPos()-R).squared_magnitude(); - float w = mexp(dpos*gamma_s); - modelVal += w*(bw+m_ParticleChemicalPotential); - w = mexp(dpos*gamma_reg_s); + float pos_w = exp(-dpos*m_SigmaInv); + + modelVal += pos_w*I_0; } neighbour = m_ParticleGrid->GetNextNeighbor(); } + modelVal += CalcI0(1.0)+m_ParticleChemicalPotential; - float energy = 2*(odfVal/m_ParticleWeight-modelVal) - (mbesseli0(1.0)+m_ParticleChemicalPotential); + float energy = 2*odfVal/m_ParticleWeight - modelVal; 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 = nullptr; int ep2 = 0; 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->GetDir(),p2->GetDir()))*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->GetPos() + (p1->GetDir() * (m_ParticleLength * ep1)); vnl_vector_fixed endPoint2 = p2->GetPos() + (p2->GetDir() * (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->GetPos() + p1->GetPos()); 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; }