diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp index 27d0cf8f76..4e95f8a1d7 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp @@ -1,443 +1,457 @@ /*=================================================================== 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 using namespace mitk; -EnergyComputer::EnergyComputer(MTRand* rgen, ItkQBallImgType* qballImage, SphereInterpolator *sp, ParticleGrid *particleGrid, float *mask, vnl_matrix_fixed rotMatrix) +EnergyComputer::EnergyComputer(ItkQBallImgType* qballImage, ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, MTRand* randGen) + : m_UseTrilinearInterpolation(true) { - m_UseTrilinearInterpolation = true; m_ParticleGrid = particleGrid; - m_RandGen = rgen; - m_RotationMatrix = rotMatrix; - m_ImageData = qballImage; - m_SphereInterpolator = sp; - m_MaskImageData = mask; + m_RandGen = randGen; + m_Image = qballImage; + m_SphereInterpolator = interpolator; + m_Mask = mask; - m_Size[0] = m_ImageData->GetLargestPossibleRegion().GetSize()[0]; - m_Size[1] = m_ImageData->GetLargestPossibleRegion().GetSize()[1]; - m_Size[2] = m_ImageData->GetLargestPossibleRegion().GetSize()[2]; + 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]; if (m_Size[0]<3 || m_Size[1]<3 || m_Size[2]<3) m_UseTrilinearInterpolation = false; - m_Spacing[0] = m_ImageData->GetSpacing()[0]; - m_Spacing[1] = m_ImageData->GetSpacing()[1]; - m_Spacing[2] = m_ImageData->GetSpacing()[2]; - - nip = QBALL_ODFSIZE; - - if (nip != sp->nverts) + m_Spacing[0] = m_Image->GetSpacing()[0]; + m_Spacing[1] = m_Image->GetSpacing()[1]; + m_Spacing[2] = m_Image->GetSpacing()[2]; + + // calculate rotation matrix + vnl_matrix temp = m_Image->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]; - cumulspatprob.resize(totsz, 0.0); - spatidx.resize(totsz, 0); + m_CumulatedSpatialProbability.resize(totsz, 0.0); + m_ActiveIndices.resize(totsz, 0); m_NumActiveVoxels = 0; - cumulspatprob[0] = 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]; - if (m_MaskImageData[idx] > 0.5) + ItkFloatImageType::IndexType index; + index[0] = x; index[1] = y; index[2] = z; + if (m_Mask->GetPixel(index) > 0.5) { - cumulspatprob[m_NumActiveVoxels+1] = cumulspatprob[m_NumActiveVoxels] + m_MaskImageData[idx]; - spatidx[m_NumActiveVoxels] = idx; + 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++) - cumulspatprob[k] /= cumulspatprob[m_NumActiveVoxels]; + m_CumulatedSpatialProbability[k] /= m_CumulatedSpatialProbability[m_NumActiveVoxels]; - fprintf(stderr,"EnergyComputer: %i active voxels found\n",m_NumActiveVoxels); + std::cout << "EnergyComputer: " << m_NumActiveVoxels << " active voxels found" << std::endl; } -void EnergyComputer::setParameters(float pwei,float pwid,float chempot_connection, float length,float curv_hardthres, float inex_balance, float chempot2, float meanv) +void EnergyComputer::SetParameters(float particleWeight, float particleWidth, float connectionPotential, float curvThres, float inexBalance, float particlePotential) { - this->chempot2 = chempot2; - meanval_sq = meanv; + m_ParticleChemicalPotential = particlePotential; + m_ConnectionPotential = connectionPotential; + m_ParticleWeight = particleWeight; - eigencon_energy = chempot_connection; - eigen_energy = 0; - particle_weight = pwei; + float bal = 1/(1+exp(-inexBalance)); + m_ExtStrength = 2*bal; + m_IntStrength = 2*(1-bal)/m_SquaredParticleLength; - float bal = 1/(1+exp(-inex_balance)); - ex_strength = 2*bal; // cleanup (todo) - in_strength = 2*(1-bal)/length/length; // cleanup (todo) - // in_strength = 0.64/length/length; // cleanup (todo) + m_CurvatureThreshold = curvThres; - particle_length_sq = length*length; - curv_hard = curv_hardthres; - - float sigma_s = pwid; + float sigma_s = particleWidth; gamma_s = 1/(sigma_s*sigma_s); - gamma_reg_s =1/(length*length/4); + gamma_reg_s =1/(m_SquaredParticleLength/4); } -void EnergyComputer::drawSpatPosition(vnl_vector_fixed& R) +// draw random position from active voxels +void EnergyComputer::DrawRandomPosition(vnl_vector_fixed& R) { float r = m_RandGen->frand(); int j; int rl = 1; int rh = m_NumActiveVoxels; while(rh != rl) { j = rl + (rh-rl)/2; - if (r < cumulspatprob[j]) + if (r < m_CumulatedSpatialProbability[j]) { rh = j; continue; } - if (r > cumulspatprob[j]) + if (r > m_CumulatedSpatialProbability[j]) { rl = j+1; continue; } break; } - R[0] = m_Spacing[0]*((float)(spatidx[rh-1] % m_Size[0]) + m_RandGen->frand()); - R[1] = m_Spacing[1]*((float)((spatidx[rh-1]/m_Size[0]) % m_Size[1]) + m_RandGen->frand()); - R[2] = m_Spacing[2]*((float)(spatidx[rh-1]/(m_Size[0]*m_Size[1])) + m_RandGen->frand()); + R[0] = m_Spacing[0]*((float)(m_ActiveIndices[rh-1] % m_Size[0]) + m_RandGen->frand()); + R[1] = m_Spacing[1]*((float)((m_ActiveIndices[rh-1]/m_Size[0]) % m_Size[1]) + m_RandGen->frand()); + R[2] = m_Spacing[2]*((float)(m_ActiveIndices[rh-1]/(m_Size[0]*m_Size[1])) + m_RandGen->frand()); } -float EnergyComputer::SpatProb(vnl_vector_fixed R) +// return spatial probability of position +float EnergyComputer::SpatProb(vnl_vector_fixed pos) { - int rx = floor(R[0]/m_Spacing[0]); - int ry = floor(R[1]/m_Spacing[1]); - int rz = floor(R[2]/m_Spacing[2]); - if (rx >= 0 && rx < m_Size[0] && ry >= 0 && ry < m_Size[1] && rz >= 0 && rz < m_Size[2]){ - return m_MaskImageData[rx + m_Size[0]* (ry + m_Size[1]*rz)]; - } + 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 &R, vnl_vector_fixed &N, float &len) +float EnergyComputer::EvaluateOdf(vnl_vector_fixed& pos, vnl_vector_fixed dir) { - const int CU = 10; - vnl_vector_fixed Rs; - float Dn = 0; - int xint,yint,zint; - - vnl_vector_fixed n; - n[0] = N[0]; - n[1] = N[1]; - n[2] = N[2]; - n = m_RotationMatrix*n; - m_SphereInterpolator->getInterpolation(n); - - for (int i=-CU; i <= CU;i++) + 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++) { - Rs = R + (N * len) * ((float)i/CU); + samplePos = pos + (dir * m_ParticleLength) * ((float)i/sampleSteps); - if (!m_UseTrilinearInterpolation) + if (!m_UseTrilinearInterpolation) // image has not enough slices to use trilinear interpolation { ItkQBallImgType::IndexType index; - index[0] = floor(R[0]/m_Spacing[0]); - index[1] = floor(R[1]/m_Spacing[1]); - index[2] = floor(R[2]/m_Spacing[2]); - if (index[0]>=0 && index[0]=0 && index[1]=0 && index[2]GetLargestPossibleRegion().IsInside(index)) { - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2]); + 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 + else // use trilinear interpolation { - float Rx = Rs[0]/m_Spacing[0]-0.5; - float Ry = Rs[1]/m_Spacing[1]-0.5; - float Rz = Rs[2]/m_Spacing[2]-0.5; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; - Dn += (m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + - m_ImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + 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; } } } - Dn *= 1.0/(2*CU+1); - return Dn; + result /= (2*sampleSteps+1); // average result over taken samples + return result; } -float EnergyComputer::computeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len, Particle *dp) +float EnergyComputer::ComputeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len, Particle *dp) { float m = SpatProb(R); if (m == 0) return -INFINITY; - float Dn = evaluateODF(R,N,len); + float Dn = EvaluateOdf(R, N); float Sn = 0; float Pn = 0; m_ParticleGrid->ComputeNeighbors(R); for (;;) { Particle *p = m_ParticleGrid->GetNextNeighbor(); if (p == 0) break; if (dp != p) { float dot = fabs(dot_product(N,p->N)); float bw = mbesseli0(dot); float dpos = (p->R-R).squared_magnitude(); float w = mexp(dpos*gamma_s); - Sn += w*(bw+chempot2); + Sn += w*(bw+m_ParticleChemicalPotential); w = mexp(dpos*gamma_reg_s); Pn += w*bw; } } - float energy = 0; - energy += 2*(Dn/particle_weight-Sn) - (mbesseli0(1.0)+chempot2); - - return energy*ex_strength; + float energy = 2*(Dn/m_ParticleWeight-Sn) - (mbesseli0(1.0)+m_ParticleChemicalPotential); + return energy*m_ExtStrength; } -float EnergyComputer::computeInternalEnergy(Particle *dp) +float EnergyComputer::ComputeInternalEnergy(Particle *dp) { - float energy = eigen_energy; + float energy = 0; if (dp->pID != -1) - energy += computeInternalEnergyConnection(dp,+1); + energy += ComputeInternalEnergyConnection(dp,+1); if (dp->mID != -1) - energy += computeInternalEnergyConnection(dp,-1); + energy += ComputeInternalEnergyConnection(dp,-1); return energy; } -float EnergyComputer::computeInternalEnergyConnection(Particle *p1,int ep1) +float EnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1) { Particle *p2 = 0; int ep2; if (ep1 == 1) p2 = m_ParticleGrid->GetParticle(p1->pID); else p2 = m_ParticleGrid->GetParticle(p1->mID); if (p2->mID == p1->ID) ep2 = -1; else if (p2->pID == p1->ID) ep2 = 1; else fprintf(stderr,"EnergyComputer_connec: Connections are inconsistent!\n"); if (p2 == 0) fprintf(stderr,"bug2"); - return computeInternalEnergyConnection(p1,ep1,p2,ep2); + return ComputeInternalEnergyConnection(p1,ep1,p2,ep2); } -float EnergyComputer::computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) +float EnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) { - if ((dot_product(p1->N,p2->N))*ep1*ep2 > -curv_hard) + if ((dot_product(p1->N,p2->N))*ep1*ep2 > -m_CurvatureThreshold) return -INFINITY; - vnl_vector_fixed R1 = p1->R + (p1->N * (p1->len * ep1)); - vnl_vector_fixed R2 = p2->R + (p2->N * (p2->len * ep2)); + vnl_vector_fixed R1 = p1->R + (p1->N * (m_ParticleLength * ep1)); + vnl_vector_fixed R2 = p2->R + (p2->N * (m_ParticleLength * ep2)); - if ((R1-R2).squared_magnitude() > particle_length_sq) + if ((R1-R2).squared_magnitude() > m_SquaredParticleLength) return -INFINITY; vnl_vector_fixed R = (p2->R + p1->R); R *= 0.5; if (SpatProb(R) == 0) return -INFINITY; float norm1 = (R1-R).squared_magnitude(); float norm2 = (R2-R).squared_magnitude(); - float energy = (eigencon_energy-norm1-norm2)*in_strength; + 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); } //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 6b73d389c9..dfee8970d4 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h @@ -1,106 +1,85 @@ /*=================================================================== 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::Image ItkFloatImageType; - vnl_matrix_fixed m_RotationMatrix; - ItkQBallImgType* m_ImageData; - vnl_vector_fixed m_Size; - vnl_vector_fixed m_Spacing; - SphereInterpolator* m_SphereInterpolator; - ParticleGrid* m_ParticleGrid; - - std::vector< float > cumulspatprob; - std::vector< int > spatidx; - - bool m_UseTrilinearInterpolation; - - float *m_MaskImageData; - int m_NumActiveVoxels; - - int w,h,d; - float voxsize_w; - float voxsize_h; - float voxsize_d; - - int w_sp,h_sp,d_sp; - float voxsize_sp_w; - float voxsize_sp_h; - float voxsize_sp_d; - MTRand* m_RandGen; - - float eigen_energy; - int nip; // number of data vertices on sphere - - - float eigencon_energy; + EnergyComputer(ItkQBallImgType* qballImage, ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, MTRand* randGen); - float chempot2; - float meanval_sq; + void SetParameters(float particleWeight, float particleWidth, float connectionPotential, float curvThres, float inexBalance, float particlePotential); - float gamma_s; - float gamma_reg_s; + void DrawRandomPosition(vnl_vector_fixed& R); - float particle_weight; - float ex_strength; - float in_strength; + float ComputeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len, Particle *dp); - float particle_length_sq; - float curv_hard; + float ComputeInternalEnergyConnection(Particle *p1,int ep1); + float ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2); - EnergyComputer(MTRand* rgen, ItkQBallImgType* qballImage, SphereInterpolator *sp, ParticleGrid *pcon, float *mask, vnl_matrix_fixed rotMatrix); + float ComputeInternalEnergy(Particle *dp); - void setParameters(float pwei,float pwid,float chempot_connection, float length,float curv_hardthres, float inex_balance, float chempot2, float meanv); +protected: - void drawSpatPosition(vnl_vector_fixed& R); - - float SpatProb(vnl_vector_fixed R); - - float evaluateODF(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len); - - float computeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len, Particle *dp); - - float computeInternalEnergy(Particle *dp); - - float computeInternalEnergyConnection(Particle *p1,int ep1); - - float computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2); + vnl_matrix_fixed m_RotationMatrix; + SphereInterpolator* m_SphereInterpolator; + ParticleGrid* m_ParticleGrid; + MTRand* m_RandGen; + 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; + + bool m_UseTrilinearInterpolation; // is deactivated if less than 3 image slices are available + + int m_NumActiveVoxels; + float m_ConnectionPotential; + float m_ParticleChemicalPotential; + float gamma_s; + float gamma_reg_s; + float m_ParticleWeight; + float m_ExtStrength; + float m_IntStrength; + float m_ParticleLength; + float m_SquaredParticleLength; + float m_CurvatureThreshold; + + 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/mitkMetropolisHastingsSampler.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp index 459633b559..48f0be2349 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp @@ -1,453 +1,450 @@ /*=================================================================== 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 "mitkMetropolisHastingsSampler.h" using namespace mitk; MetropolisHastingsSampler::MetropolisHastingsSampler(ParticleGrid* grid, EnergyComputer* enComp, MTRand* randGen, float curvThres) : m_AcceptedProposals(0) , m_ExTemp(0.01) , m_BirthProb(0.25) , m_DeathProb(0.05) , m_ShiftProb(0.15) , m_OptShiftProb(0.1) , m_ConnectionProb(0.45) , m_TractProb(0.5) , m_DelProb(0.1) , m_ChempotParticle(0.0) { m_RandGen = randGen; m_ParticleGrid = grid; m_EnergyComputer = enComp; m_ParticleLength = m_ParticleGrid->m_ParticleLength; m_DistanceThreshold = m_ParticleLength*m_ParticleLength; m_Sigma = m_ParticleLength/8.0; m_Gamma = 1/(m_Sigma*m_Sigma*2); m_Z = pow(2*M_PI*m_Sigma,3.0/2.0)*(M_PI*m_Sigma/m_ParticleLength); m_CurvatureThreshold = curvThres; // ??? m_StopProb = exp(-1/m_TractProb); } void MetropolisHastingsSampler::SetTemperature(float val) { m_InTemp = val; m_Density = exp(-m_ChempotParticle/m_InTemp); } vnl_vector_fixed MetropolisHastingsSampler::DistortVector(float sigma, vnl_vector_fixed& vec) { vec[0] += sigma*m_RandGen->frandn(); vec[1] += sigma*m_RandGen->frandn(); vec[2] += sigma*m_RandGen->frandn(); } vnl_vector_fixed MetropolisHastingsSampler::GetRandomDirection() { vnl_vector_fixed vec; vec[0] += m_RandGen->frandn(); vec[1] += m_RandGen->frandn(); vec[2] += m_RandGen->frandn(); vec.normalize(); return vec; } void MetropolisHastingsSampler::MakeProposal() { float randnum = m_RandGen->frand(); // Birth Proposal if (randnum < m_BirthProb) { vnl_vector_fixed R; - m_EnergyComputer->drawSpatPosition(R); + m_EnergyComputer->DrawRandomPosition(R); vnl_vector_fixed N = GetRandomDirection(); - float len = m_ParticleLength; Particle prop; prop.R = R; prop.N = N; - prop.len = len; float prob = m_Density * m_DeathProb /((m_BirthProb)*(m_ParticleGrid->m_NumParticles+1)); - float ex_energy = m_EnergyComputer->computeExternalEnergy(R,N,len,0); - float in_energy = m_EnergyComputer->computeInternalEnergy(&prop); + float ex_energy = m_EnergyComputer->ComputeExternalEnergy(R,N,m_ParticleLength,0); + float in_energy = m_EnergyComputer->ComputeInternalEnergy(&prop); prob *= exp((in_energy/m_InTemp+ex_energy/m_ExTemp)) ; if (prob > 1 || m_RandGen->frand() < prob) { Particle *p = m_ParticleGrid->NewParticle(R); if (p!=0) { p->R = R; p->N = N; - p->len = len; m_AcceptedProposals++; } } } // Death Proposal else if (randnum < m_BirthProb+m_DeathProb) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *dp = m_ParticleGrid->GetParticle(pnum); if (dp->pID == -1 && dp->mID == -1) { - float ex_energy = m_EnergyComputer->computeExternalEnergy(dp->R,dp->N,dp->len,dp); - float in_energy = m_EnergyComputer->computeInternalEnergy(dp); + float ex_energy = m_EnergyComputer->ComputeExternalEnergy(dp->R,dp->N,m_ParticleLength,dp); + float in_energy = m_EnergyComputer->ComputeInternalEnergy(dp); float prob = m_ParticleGrid->m_NumParticles * (m_BirthProb) /(m_Density*m_DeathProb); //*SpatProb(dp->R); prob *= exp(-(in_energy/m_InTemp+ex_energy/m_ExTemp)) ; if (prob > 1 || m_RandGen->frand() < prob) { m_ParticleGrid->RemoveParticle(pnum); m_AcceptedProposals++; } } } } // Shift Proposal else if (randnum < m_BirthProb+m_DeathProb+m_ShiftProb) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = m_ParticleGrid->GetParticle(pnum); Particle prop_p = *p; DistortVector(m_Sigma, prop_p.R); - DistortVector(m_Sigma/(2*p->len), prop_p.N); + DistortVector(m_Sigma/(2*m_ParticleLength), prop_p.N); prop_p.N.normalize(); - float ex_energy = m_EnergyComputer->computeExternalEnergy(prop_p.R,prop_p.N,p->len,p) - - m_EnergyComputer->computeExternalEnergy(p->R,p->N,p->len,p); - float in_energy = m_EnergyComputer->computeInternalEnergy(&prop_p) - m_EnergyComputer->computeInternalEnergy(p); + float ex_energy = m_EnergyComputer->ComputeExternalEnergy(prop_p.R,prop_p.N,m_ParticleLength,p) + - m_EnergyComputer->ComputeExternalEnergy(p->R,p->N,m_ParticleLength,p); + float in_energy = m_EnergyComputer->ComputeInternalEnergy(&prop_p) - m_EnergyComputer->ComputeInternalEnergy(p); float prob = exp(ex_energy/m_ExTemp+in_energy/m_InTemp); if (m_RandGen->frand() < prob) { vnl_vector_fixed Rtmp = p->R; vnl_vector_fixed Ntmp = p->N; p->R = prop_p.R; p->N = prop_p.N; if (!m_ParticleGrid->TryUpdateGrid(pnum)) { p->R = Rtmp; p->N = Ntmp; } m_AcceptedProposals++; } } } // Optimal Shift Proposal else if (randnum < m_BirthProb+m_DeathProb+m_ShiftProb+m_OptShiftProb) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = m_ParticleGrid->GetParticle(pnum); bool no_proposal = false; Particle prop_p = *p; if (p->pID != -1 && p->mID != -1) { Particle *plus = m_ParticleGrid->GetParticle(p->pID); int ep_plus = (plus->pID == p->ID)? 1 : -1; Particle *minus = m_ParticleGrid->GetParticle(p->mID); int ep_minus = (minus->pID == p->ID)? 1 : -1; - prop_p.R = (plus->R + plus->N * (plus->len * ep_plus) + minus->R + minus->N * (minus->len * ep_minus)); + prop_p.R = (plus->R + plus->N * (m_ParticleLength * ep_plus) + minus->R + minus->N * (m_ParticleLength * ep_minus)); prop_p.R *= 0.5; prop_p.N = plus->R - minus->R; prop_p.N.normalize(); } else if (p->pID != -1) { Particle *plus = m_ParticleGrid->GetParticle(p->pID); int ep_plus = (plus->pID == p->ID)? 1 : -1; - prop_p.R = plus->R + plus->N * (plus->len * ep_plus * 2); + prop_p.R = plus->R + plus->N * (m_ParticleLength * ep_plus * 2); prop_p.N = plus->N; } else if (p->mID != -1) { Particle *minus = m_ParticleGrid->GetParticle(p->mID); int ep_minus = (minus->pID == p->ID)? 1 : -1; - prop_p.R = minus->R + minus->N * (minus->len * ep_minus * 2); + prop_p.R = minus->R + minus->N * (m_ParticleLength * ep_minus * 2); prop_p.N = minus->N; } else no_proposal = true; if (!no_proposal) { float cos = dot_product(prop_p.N, p->N); float p_rev = exp(-((prop_p.R-p->R).squared_magnitude() + (1-cos*cos))*m_Gamma)/m_Z; - float ex_energy = m_EnergyComputer->computeExternalEnergy(prop_p.R,prop_p.N,p->len,p) - - m_EnergyComputer->computeExternalEnergy(p->R,p->N,p->len,p); - float in_energy = m_EnergyComputer->computeInternalEnergy(&prop_p) - m_EnergyComputer->computeInternalEnergy(p); + float ex_energy = m_EnergyComputer->ComputeExternalEnergy(prop_p.R,prop_p.N,m_ParticleLength,p) + - m_EnergyComputer->ComputeExternalEnergy(p->R,p->N,m_ParticleLength,p); + float in_energy = m_EnergyComputer->ComputeInternalEnergy(&prop_p) - m_EnergyComputer->ComputeInternalEnergy(p); float prob = exp(ex_energy/m_ExTemp+in_energy/m_InTemp)*m_ShiftProb*p_rev/(m_OptShiftProb+m_ShiftProb*p_rev); if (m_RandGen->frand() < prob) { vnl_vector_fixed Rtmp = p->R; vnl_vector_fixed Ntmp = p->N; p->R = prop_p.R; p->N = prop_p.N; if (!m_ParticleGrid->TryUpdateGrid(pnum)) { p->R = Rtmp; p->N = Ntmp; } m_AcceptedProposals++; } } } } else { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = m_ParticleGrid->GetParticle(pnum); EndPoint P; P.p = p; P.ep = (m_RandGen->frand() > 0.5)? 1 : -1; RemoveAndSaveTrack(P); if (m_BackupTrack.m_Probability != 0) { MakeTrackProposal(P); float prob = (m_ProposalTrack.m_Energy-m_BackupTrack.m_Energy)/m_InTemp ; prob = exp(prob)*(m_BackupTrack.m_Probability * pow(m_DelProb,m_ProposalTrack.m_Length)) /(m_ProposalTrack.m_Probability * pow(m_DelProb,m_BackupTrack.m_Length)); if (m_RandGen->frand() < prob) { ImplementTrack(m_ProposalTrack); m_AcceptedProposals++; } else { ImplementTrack(m_BackupTrack); } } else ImplementTrack(m_BackupTrack); } } } void MetropolisHastingsSampler::ImplementTrack(Track &T) { for (int k = 1; k < T.m_Length;k++) m_ParticleGrid->CreateConnection(T.track[k-1].p,T.track[k-1].ep,T.track[k].p,-T.track[k].ep); } void MetropolisHastingsSampler::RemoveAndSaveTrack(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; m_BackupTrack.track[cnt] = Current; EndPoint Next; for (;;) { Next.p = 0; if (Current.ep == 1) { if (Current.p->pID != -1) { Next.p = m_ParticleGrid->GetParticle(Current.p->pID); Current.p->pID = -1; m_ParticleGrid->m_NumConnections--; } } else if (Current.ep == -1) { if (Current.p->mID != -1) { Next.p = m_ParticleGrid->GetParticle(Current.p->mID); Current.p->mID = -1; m_ParticleGrid->m_NumConnections--; } } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 3\n"); break; } if (Next.p == 0) // no successor { Next.ep = 0; // mark as empty successor break; } else { if (Next.p->pID == Current.p->ID) { Next.p->pID = -1; Next.ep = 1; } else if (Next.p->mID == Current.p->ID) { Next.p->mID = -1; Next.ep = -1; } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 4\n"); break; } } ComputeEndPointProposalDistribution(Current); AccumProb *= (m_SimpSamp.probFor(Next)); if (Next.p == 0) // no successor -> break break; - energy += m_EnergyComputer->computeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); + energy += m_EnergyComputer->ComputeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); Current = Next; Current.ep *= -1; cnt++; m_BackupTrack.track[cnt] = Current; if (m_RandGen->rand() > m_DelProb) break; } m_BackupTrack.m_Energy = energy; m_BackupTrack.m_Probability = AccumProb; m_BackupTrack.m_Length = cnt+1; } void MetropolisHastingsSampler::MakeTrackProposal(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; m_ProposalTrack.track[cnt++] = Current; Current.p->label = 1; for (;;) { // next candidate is already connected if ((Current.ep == 1 && Current.p->pID != -1) || (Current.ep == -1 && Current.p->mID != -1)) break; // track too long // if (cnt > 250) // break; ComputeEndPointProposalDistribution(Current); int k = m_SimpSamp.draw(m_RandGen->frand()); // stop tracking proposed if (k==0) break; EndPoint Next = m_SimpSamp.objs[k]; float probability = m_SimpSamp.probFor(k); // accumulate energy and proposal distribution - energy += m_EnergyComputer->computeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); + energy += m_EnergyComputer->ComputeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); AccumProb *= probability; // track to next endpoint Current = Next; Current.ep *= -1; Current.p->label = 1; // put label to avoid loops m_ProposalTrack.track[cnt++] = Current; } m_ProposalTrack.m_Energy = energy; m_ProposalTrack.m_Probability = AccumProb; m_ProposalTrack.m_Length = cnt; // clear labels for (int j = 0; j < m_ProposalTrack.m_Length;j++) m_ProposalTrack.track[j].p->label = 0; } void MetropolisHastingsSampler::ComputeEndPointProposalDistribution(EndPoint P) { Particle *p = P.p; int ep = P.ep; float dist,dot; - vnl_vector_fixed R = p->R + (p->N * (ep*p->len) ); + vnl_vector_fixed R = p->R + (p->N * (ep*m_ParticleLength) ); m_ParticleGrid->ComputeNeighbors(R); m_SimpSamp.clear(); m_SimpSamp.add(m_StopProb,EndPoint(0,0)); for (;;) { Particle *p2 = m_ParticleGrid->GetNextNeighbor(); if (p2 == 0) break; if (p!=p2 && p2->label == 0) { if (p2->mID == -1) { - dist = (p2->R - p2->N * p2->len - R).squared_magnitude(); + dist = (p2->R - p2->N * m_ParticleLength - R).squared_magnitude(); if (dist < m_DistanceThreshold) { dot = dot_product(p2->N,p->N) * ep; if (dot > m_CurvatureThreshold) { - float en = m_EnergyComputer->computeInternalEnergyConnection(p,ep,p2,-1); + float en = m_EnergyComputer->ComputeInternalEnergyConnection(p,ep,p2,-1); m_SimpSamp.add(exp(en/m_TractProb),EndPoint(p2,-1)); } } } if (p2->pID == -1) { - dist = (p2->R + p2->N * p2->len - R).squared_magnitude(); + dist = (p2->R + p2->N * m_ParticleLength - R).squared_magnitude(); if (dist < m_DistanceThreshold) { dot = dot_product(p2->N,p->N) * (-ep); if (dot > m_CurvatureThreshold) { - float en = m_EnergyComputer->computeInternalEnergyConnection(p,ep,p2,+1); + float en = m_EnergyComputer->ComputeInternalEnergyConnection(p,ep,p2,+1); m_SimpSamp.add(exp(en/m_TractProb),EndPoint(p2,+1)); } } } } } } int MetropolisHastingsSampler::GetNumAcceptedProposals() { return m_AcceptedProposals; } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h index f1741906a4..5dbdf1311c 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h @@ -1,75 +1,74 @@ /*=================================================================== 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 _PARTICLE #define _PARTICLE #include #include namespace mitk { class MitkDiffusionImaging_EXPORT Particle { public: Particle() { label = 0; pID = -1; mID = -1; } ~Particle() { } vnl_vector_fixed R; vnl_vector_fixed N; - float len; int gridindex; // index in the grid where it is living int ID; int pID; int mID; unsigned char label; }; class MitkDiffusionImaging_EXPORT EndPoint { public: EndPoint() {} EndPoint(Particle *p,int ep) { this->p = p; this->ep = ep; } Particle *p; int ep; inline bool operator==(EndPoint P) { return (P.p == p) && (P.ep == ep); } }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index 78e00d18f3..9eba2dde3a 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,389 +1,373 @@ /*=================================================================== 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 // ITK #include -#pragma GCC visibility push(default) -#include -#pragma GCC visibility pop +#include // MISC -#include #include #include namespace itk{ template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::GibbsTrackingFilter(): - m_TempStart(0.1), - m_TempEnd(0.001), - m_NumIt(500000), + m_StartTemperature(0.1), + m_EndTemperature(0.001), + m_Iterations(500000), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), - m_ChempotConnection(10), + m_ConnectionPotential(10), m_InexBalance(0), - m_Chempot2(0.2), - m_FiberLength(10), + 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_Memory(0), m_ProposalAcceptance(0), - m_CurvatureHardThreshold(0.7), - m_Meanval_sq(0.0), + m_CurvatureThreshold(0.7), m_DuplicateImage(true) { } 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 > bool GibbsTrackingFilter< ItkQBallImageType > ::EstimateParticleWeight() { MITK_INFO << "itkGibbsTrackingFilter: estimating particle weight"; typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(m_QBallImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); ItkFloatImageType::Pointer gfaImage = gfaFilter->GetOutput(); float samplingStart = 1.0; float samplingStop = 0.66; // GFA iterator typedef ImageRegionIterator< ItkFloatImageType > GfaIteratorType; GfaIteratorType gfaIt(gfaImage, gfaImage->GetLargestPossibleRegion() ); // Mask iterator typedef ImageRegionConstIterator< ItkFloatImageType > MaskIteratorType; MaskIteratorType mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); // Input iterator typedef ImageRegionConstIterator< ItkQBallImageType > InputIteratorType; InputIteratorType it(m_QBallImage, m_QBallImage->GetLargestPossibleRegion() ); float upper = 0; int count = 0; for(float thr=samplingStart; thr>samplingStop; thr-=0.01) { it.GoToBegin(); mit.GoToBegin(); gfaIt.GoToBegin(); while( !gfaIt.IsAtEnd() ) { if(gfaIt.Get()>thr && mit.Get()>0) { itk::OrientationDistributionFunction odf(it.Get().GetDataPointer()); upper += odf.GetMaxValue()-odf.GetMeanValue(); ++count; } ++it; ++mit; ++gfaIt; } } if (count>0) upper /= count; else return false; m_ParticleWeight = upper/6; return true; } // 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(); } - // image sizes and spacing - int imgSize[4] = { QBALL_ODFSIZE, - m_QBallImage->GetLargestPossibleRegion().GetSize().GetElement(0), - m_QBallImage->GetLargestPossibleRegion().GetSize().GetElement(1), - m_QBallImage->GetLargestPossibleRegion().GetSize().GetElement(2)}; - double spacing[3] = {m_QBallImage->GetSpacing().GetElement(0),m_QBallImage->GetSpacing().GetElement(1),m_QBallImage->GetSpacing().GetElement(2)}; - - // calculate rotation matrix - vnl_matrix temp = m_QBallImage->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)){ - MITK_INFO << "itkGibbsTrackingFilter: image direction is not a rotation matrix. Tracking not possible!"; - m_AbortTracking = true; - } - - // generate local working copy of QBall image - typename ItkQBallImageType::Pointer qballImage; + // generate local working copy of QBall image (if not disabled) if (m_DuplicateImage) { typedef itk::ImageDuplicator< ItkQBallImageType > DuplicateFilterType; typename DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); duplicator->SetInputImage( m_QBallImage ); duplicator->Update(); - qballImage = duplicator->GetOutput(); - } - else - { - qballImage = m_QBallImage; + m_QBallImage = duplicator->GetOutput(); } // perform mean subtraction on odfs typedef ImageRegionIterator< ItkQBallImageType > InputIteratorType; - InputIteratorType it(qballImage, qballImage->GetLargestPossibleRegion() ); + 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; } - // mask image - int maskImageSize[3]; - float *mask; - if(m_MaskImage.IsNull()) - { - m_MaskImage = ItkFloatImageType::New(); - m_MaskImage->SetSpacing( qballImage->GetSpacing() ); // Set the image spacing - m_MaskImage->SetOrigin( qballImage->GetOrigin() ); // Set the image origin - m_MaskImage->SetDirection( qballImage->GetDirection() ); // Set the image direction - m_MaskImage->SetRegions( qballImage->GetLargestPossibleRegion()); - m_MaskImage->Allocate(); - m_MaskImage->FillBuffer(1.0); - } - mask = (float*) m_MaskImage->GetBufferPointer(); - maskImageSize[0] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(0); - maskImageSize[1] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(1); - maskImageSize[2] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(2); + // check if mask image is given if it needs resampling + PrepareMaskImage(); - // get paramters + // prepare parameters float minSpacing; - if(spacing[0]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 = spacing[2]; + 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) if (!EstimateParticleWeight()) { - MITK_INFO << "itkGibbsTrackingFilter: could not estimate particle weight!"; + MITK_INFO << "itkGibbsTrackingFilter: could not estimate particle weight. using default value."; m_ParticleWeight = 0.0001; } - m_CurrentStep = 0; - m_Memory = 0; - - float alpha = log(m_TempEnd/m_TempStart); - m_Steps = m_NumIt/10000; + float alpha = log(m_EndTemperature/m_StartTemperature); + m_Steps = m_Iterations/10000; if (m_Steps<10) m_Steps = 10; - if (m_Steps>m_NumIt) + if (m_Steps>m_Iterations) { MITK_INFO << "itkGibbsTrackingFilter: 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)); - if (m_CurvatureHardThreshold < mitk::eps) - m_CurvatureHardThreshold = 0; - unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); - + // seed random generators MTRand randGen(1); srand(1); + // load sphere interpolator to evaluate the ODFs SphereInterpolator* interpolator = LoadSphereInterpolator(); - MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; + // initialize the actual tracking components (ParticleGrid, Metropolis Hastings Sampler and Energy Computer) ParticleGrid* particleGrid = new ParticleGrid(m_MaskImage, m_ParticleLength); - EnergyComputer encomp(&randGen, qballImage, interpolator, particleGrid, mask, directionMatrix); - encomp.setParameters(m_ParticleWeight,m_ParticleWidth,m_ChempotConnection*m_ParticleLength*m_ParticleLength,m_ParticleLength,m_CurvatureHardThreshold,m_InexBalance,m_Chempot2, m_Meanval_sq); + EnergyComputer encomp(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_CurvatureHardThreshold); + MetropolisHastingsSampler sampler(particleGrid, &encomp, &randGen, m_CurvatureThreshold); // main loop m_NumAcceptedFibers = 0; for( m_CurrentStep = 1; m_CurrentStep <= m_Steps; m_CurrentStep++ ) { - m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/m_NumIt; - m_NumParticles = particleGrid->m_NumParticles; - m_NumConnections = particleGrid->m_NumConnections; - - MITK_INFO << "itkGibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; - MITK_INFO << "itkGibbsTrackingFilter: particles: " << m_NumParticles; - MITK_INFO << "itkGibbsTrackingFilter: connections: " << m_NumConnections; - MITK_INFO << "itkGibbsTrackingFilter: progress: " << 100*(float)m_CurrentStep/m_Steps << "%"; - - float temperature = m_TempStart * exp(alpha*(((1.0)*m_CurrentStep)/((1.0)*m_Steps))); - sampler->SetTemperature(temperature); + // 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(); + sampler.MakeProposal(); - if (m_BuildFibers) + if (m_BuildFibers || (i==singleIts-1 && m_CurrentStep==m_Steps)) { - m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/m_NumIt; + m_ProposalAcceptance = (float)sampler.GetNumAcceptedProposals()/m_Iterations; m_NumParticles = particleGrid->m_NumParticles; m_NumConnections = particleGrid->m_NumConnections; FiberBuilder fiberBuilder(particleGrid, m_MaskImage); - m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); + m_FiberPolyData = fiberBuilder.iterate(m_MinFiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); m_BuildFibers = false; } } + + m_ProposalAcceptance = (float)sampler.GetNumAcceptedProposals()/m_Iterations; + m_NumParticles = particleGrid->m_NumParticles; + m_NumConnections = particleGrid->m_NumConnections; + MITK_INFO << "itkGibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; + MITK_INFO << "itkGibbsTrackingFilter: particles: " << m_NumParticles; + MITK_INFO << "itkGibbsTrackingFilter: connections: " << m_NumConnections; + MITK_INFO << "itkGibbsTrackingFilter: progress: " << 100*(float)m_CurrentStep/m_Steps << "%"; } - FiberBuilder fiberBuilder(particleGrid, m_MaskImage); - m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); - m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); delete interpolator; - delete sampler; delete particleGrid; m_AbortTracking = true; m_BuildFibers = false; MITK_INFO << "itkGibbsTrackingFilter: done generate data"; } +template< class ItkQBallImageType > +void GibbsTrackingFilter< ItkQBallImageType >::PrepareMaskImage() +{ + if(m_MaskImage.IsNull()) + { + MITK_INFO << "itkGibbsTrackingFilter: 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 << "itkGibbsTrackingFilter: 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 << "itkGibbsTrackingFilter: resampling finished"; + } +} + template< class ItkQBallImageType > SphereInterpolator* GibbsTrackingFilter< ItkQBallImageType >::LoadSphereInterpolator() { QString applicationDir = QCoreApplication::applicationDirPath(); applicationDir.append("/"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); applicationDir.append("../"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); applicationDir.append("../../"); mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); std::string lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTBaryCoords.bin"); ifstream BaryCoords; BaryCoords.open(lutPath.c_str(), ios::in | ios::binary); float* coords; if (BaryCoords.is_open()) { float tmp; coords = new float [1630818]; BaryCoords.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { BaryCoords.read((char *)&tmp, sizeof(tmp)); coords[i] = tmp; } BaryCoords.close(); } else { MITK_INFO << "itkGibbsTrackingFilter: unable to open barycoords file"; m_AbortTracking = true; } ifstream Indices; lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTIndices.bin"); Indices.open(lutPath.c_str(), ios::in | ios::binary); int* ind; if (Indices.is_open()) { int tmp; ind = new int [1630818]; Indices.seekg (0, ios::beg); for (int i=0; i<1630818; i++) { Indices.read((char *)&tmp, 4); ind[i] = tmp; } Indices.close(); } else { MITK_INFO << "itkGibbsTrackingFilter: unable to open indices file"; m_AbortTracking = true; } // initialize sphere interpolator with lookuptables return new SphereInterpolator(coords, ind, QBALL_ODFSIZE, 301, 0.5); } } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h index d198de102b..60ecb41862 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h @@ -1,157 +1,135 @@ /*=================================================================== 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 itkGibbsTrackingFilter_h #define itkGibbsTrackingFilter_h // MITK #include // ITK #include #include #include // VTK #include #include #include #include #include namespace itk{ template< class ItkQBallImageType > class GibbsTrackingFilter : public ProcessObject { public: typedef GibbsTrackingFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; itkNewMacro(Self) itkTypeMacro( GibbsTrackingFilter, ProcessObject ) - typedef Image< DiffusionTensor3D, 3 > ItkTensorImage; - typedef typename ItkQBallImageType::Pointer ItkQBallImageTypePointer; - typedef Image< float, 3 > ItkFloatImageType; - typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; - - itkSetMacro( TempStart, float ) - itkGetMacro( TempStart, float ) - - itkSetMacro( TempEnd, float ) - itkGetMacro( TempEnd, float ) - - itkSetMacro( NumIt, unsigned long ) - itkGetMacro( NumIt, unsigned long ) - + typedef Image< DiffusionTensor3D, 3 > ItkTensorImage; + typedef typename ItkQBallImageType::Pointer ItkQBallImageTypePointer; + typedef Image< float, 3 > ItkFloatImageType; + typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; + + // parameter setter + itkSetMacro( StartTemperature, float ) + itkSetMacro( EndTemperature, float ) + itkSetMacro( Iterations, unsigned long ) itkSetMacro( ParticleWeight, float ) - itkGetMacro( ParticleWeight, float ) - - /** width of particle sigma (std-dev of gaussian around center) **/ itkSetMacro( ParticleWidth, float ) - itkGetMacro( ParticleWidth, float ) - - /** length of particle from midpoint to ends **/ itkSetMacro( ParticleLength, float ) - itkGetMacro( ParticleLength, float ) - - itkSetMacro( ChempotConnection, float ) - itkGetMacro( ChempotConnection, float ) - + itkSetMacro( ConnectionPotential, float ) itkSetMacro( InexBalance, float ) - itkGetMacro( InexBalance, float ) - - itkSetMacro( Chempot2, float ) - itkGetMacro( Chempot2, float ) - - itkSetMacro( FiberLength, int ) - itkGetMacro( FiberLength, int ) - + itkSetMacro( ParticlePotential, float ) + itkSetMacro( MinFiberLength, int ) itkSetMacro( AbortTracking, bool ) - itkGetMacro( AbortTracking, bool ) + itkSetMacro( CurvatureThreshold, float) + itkSetMacro( DuplicateImage, bool ) - itkSetMacro( CurrentStep, unsigned long ) + // getter + itkGetMacro( ParticleWeight, float ) + itkGetMacro( ParticleWidth, float ) + itkGetMacro( ParticleLength, float ) itkGetMacro( CurrentStep, unsigned long ) - - itkSetMacro( CurvatureHardThreshold, float) - itkGetMacro( CurvatureHardThreshold, float) - - itkGetMacro(NumParticles, unsigned long) - itkGetMacro(NumConnections, unsigned long) - itkGetMacro(NumAcceptedFibers, int) - itkGetMacro(ProposalAcceptance, float) - itkGetMacro(Steps, unsigned int) + itkGetMacro( NumParticles, int ) + itkGetMacro( NumConnections, int ) + itkGetMacro( NumAcceptedFibers, int ) + itkGetMacro( ProposalAcceptance, float ) + itkGetMacro( Steps, unsigned int) // input data itkSetMacro(QBallImage, typename ItkQBallImageType::Pointer) itkSetMacro(MaskImage, ItkFloatImageType::Pointer) itkSetMacro(TensorImage, ItkTensorImage::Pointer) void GenerateData(); virtual void Update(){ this->GenerateData(); } FiberPolyDataType GetFiberBundle(); protected: GibbsTrackingFilter(); virtual ~GibbsTrackingFilter(); bool EstimateParticleWeight(); SphereInterpolator* LoadSphereInterpolator(); + void PrepareMaskImage(); // Input Images typename ItkQBallImageType::Pointer m_QBallImage; typename ItkFloatImageType::Pointer m_MaskImage; typename ItkTensorImage::Pointer m_TensorImage; // Tracking parameters - float m_TempStart; // Start temperature - float m_TempEnd; // End temperature - unsigned long m_NumIt; // Total number of iterations - unsigned long m_CurrentStep; // current tracking step - float m_ParticleWeight; // w (unitless) - float m_ParticleWidth; //sigma (mm) - float m_ParticleLength; // ell (mm) - float m_ChempotConnection; // gross L (chemisches potential) - float m_InexBalance; // gewichtung zwischen den lambdas; -5 ... 5 -> nur intern ... nur extern,default 0 - float m_Chempot2; // typischerweise 0 - int m_FiberLength; - bool m_AbortTracking; - int m_NumAcceptedFibers; - volatile bool m_BuildFibers; - unsigned int m_Steps; - float m_Memory; - float m_ProposalAcceptance; - float m_CurvatureHardThreshold; - float m_Meanval_sq; - bool m_DuplicateImage; - - FiberPolyDataType m_FiberPolyData; - unsigned long m_NumParticles; - unsigned long m_NumConnections; + float m_StartTemperature; // Start temperature + float m_EndTemperature; // End temperature + unsigned long m_Iterations; // Total number of iterations + unsigned long m_CurrentStep; // current tracking step + float m_ParticleWeight; // w (unitless) + float m_ParticleWidth; // sigma (mm) + float m_ParticleLength; // l (mm) + float m_ConnectionPotential; // gross L (chemisches potential, default 10) + float m_InexBalance; // gewichtung zwischen den lambdas; -5 ... 5 -> nur intern ... nur extern,default 0 + float m_ParticlePotential; // default 0.2 + int m_MinFiberLength; // discard all fibers shortan than the specified length in mm + bool m_AbortTracking; // set flag to abort tracking + int m_NumAcceptedFibers; // number of reconstructed fibers generated by the FiberBuilder + volatile bool m_BuildFibers; // set flag to generate fibers from particle grid + unsigned int m_Steps; // number of temperature decrease steps + float m_ProposalAcceptance; // proposal acceptance rate (0-1) + float m_CurvatureThreshold; // curvature threshold in radians (1 -> no curvature is accepted, -1 all curvature angles are accepted) + bool m_DuplicateImage; // generates a working copy of the qball image so that the original image won't be changed by the mean subtraction + int m_NumParticles; // current number of particles in grid + int m_NumConnections; // current number of connections between particles in grid + + FiberPolyDataType m_FiberPolyData; // container for reconstructed fibers }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkGibbsTrackingFilter.cpp" #endif #endif diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp index 656c2e207d..adac6af8ac 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp @@ -1,750 +1,740 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkGibbsTrackingView.h" #include // Qt #include #include #include // MITK #include #include #include #include #include // ITK #include #include #include // MISC #include QmitkTrackingWorker::QmitkTrackingWorker(QmitkGibbsTrackingView* view) - : m_View(view) + : m_View(view) { } void QmitkTrackingWorker::run() { - m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); - -// MITK_INFO << "Resampling mask images"; -// // setup resampler -// typedef itk::ResampleImageFilter ResamplerType; -// ResamplerType::Pointer resampler = ResamplerType::New(); -// resampler->SetOutputSpacing( m_View->m_ItkQBallImage->GetSpacing() ); -// resampler->SetOutputOrigin( m_View->m_ItkQBallImage->GetOrigin() ); -// resampler->SetOutputDirection( m_View->m_ItkQBallImage->GetDirection() ); -// resampler->SetSize( m_View->m_ItkQBallImage->GetLargestPossibleRegion().GetSize() ); - -// // resample mask image -// resampler->SetInput( m_View->m_MaskImage ); -// resampler->SetDefaultPixelValue(0); -// resampler->Update(); -// m_View->m_MaskImage = resampler->GetOutput(); - - m_View->m_GlobalTracker->SetQBallImage(m_View->m_ItkQBallImage); - m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); - m_View->m_GlobalTracker->SetTempStart((float)m_View->m_Controls->m_StartTempSlider->value()/100); - m_View->m_GlobalTracker->SetTempEnd((float)m_View->m_Controls->m_EndTempSlider->value()/10000); - m_View->m_GlobalTracker->SetNumIt(m_View->m_Iterations); - m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); - m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); - m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); - m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); - m_View->m_GlobalTracker->SetFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); - m_View->m_GlobalTracker->SetCurvatureHardThreshold(cos((float)m_View->m_Controls->m_CurvatureThresholdSlider->value()*3.14159265/180)); - - m_View->m_GlobalTracker->Update(); - m_View->m_TrackingThread.quit(); + m_View->m_GlobalTracker = QmitkGibbsTrackingView::GibbsTrackingFilterType::New(); + + m_View->m_GlobalTracker->SetQBallImage(m_View->m_ItkQBallImage); + m_View->m_GlobalTracker->SetTensorImage(m_View->m_ItkTensorImage); + m_View->m_GlobalTracker->SetMaskImage(m_View->m_MaskImage); + m_View->m_GlobalTracker->SetStartTemperature((float)m_View->m_Controls->m_StartTempSlider->value()/100); + m_View->m_GlobalTracker->SetEndTemperature((float)m_View->m_Controls->m_EndTempSlider->value()/10000); + m_View->m_GlobalTracker->SetIterations(m_View->m_Iterations); + m_View->m_GlobalTracker->SetParticleWeight((float)m_View->m_Controls->m_ParticleWeightSlider->value()/10000); + m_View->m_GlobalTracker->SetParticleWidth((float)(m_View->m_Controls->m_ParticleWidthSlider->value())/10); + m_View->m_GlobalTracker->SetParticleLength((float)(m_View->m_Controls->m_ParticleLengthSlider->value())/10); + m_View->m_GlobalTracker->SetInexBalance((float)m_View->m_Controls->m_InExBalanceSlider->value()/10); + m_View->m_GlobalTracker->SetMinFiberLength(m_View->m_Controls->m_FiberLengthSlider->value()); + m_View->m_GlobalTracker->SetCurvatureThreshold(cos((float)m_View->m_Controls->m_CurvatureThresholdSlider->value()*3.14159265/180)); + + m_View->m_GlobalTracker->Update(); + m_View->m_TrackingThread.quit(); } const std::string QmitkGibbsTrackingView::VIEW_ID = -"org.mitk.views.gibbstracking"; + "org.mitk.views.gibbstracking"; QmitkGibbsTrackingView::QmitkGibbsTrackingView() - : QmitkFunctionality() - , m_Controls( 0 ) - , m_MultiWidget( NULL ) - , m_ThreadIsRunning(false) - , m_GlobalTracker(NULL) - , m_QBallImage(NULL) - , m_MaskImage(NULL) - , m_QBallImageNode(NULL) - , m_ItkQBallImage(NULL) - , m_FiberBundleNode(NULL) - , m_MaskImageNode(NULL) - , m_TrackingWorker(this) - , m_Iterations(10000000) - , m_LastStep(0) -{ - m_TrackingWorker.moveToThread(&m_TrackingThread); - connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); - connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); - connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); - connect(&m_TrackingThread, SIGNAL(terminated()), this, SLOT(AfterThread())); - m_TrackingTimer = new QTimer(this); + : QmitkFunctionality() + , m_Controls( 0 ) + , m_MultiWidget( NULL ) + , m_ThreadIsRunning(false) + , m_GlobalTracker(NULL) + , m_QBallImage(NULL) + , m_MaskImage(NULL) + , m_ImageNode(NULL) + , m_ItkQBallImage(NULL) + , m_ItkTensorImage(NULL) + , m_FiberBundleNode(NULL) + , m_MaskImageNode(NULL) + , m_TrackingWorker(this) + , m_Iterations(10000000) + , m_LastStep(0) +{ + m_TrackingWorker.moveToThread(&m_TrackingThread); + connect(&m_TrackingThread, SIGNAL(started()), this, SLOT(BeforeThread())); + connect(&m_TrackingThread, SIGNAL(started()), &m_TrackingWorker, SLOT(run())); + connect(&m_TrackingThread, SIGNAL(finished()), this, SLOT(AfterThread())); + connect(&m_TrackingThread, SIGNAL(terminated()), this, SLOT(AfterThread())); + m_TrackingTimer = new QTimer(this); } QmitkGibbsTrackingView::~QmitkGibbsTrackingView() { - delete m_TrackingTimer; + delete m_TrackingTimer; } // update tracking status and generate fiber bundle void QmitkGibbsTrackingView::TimerUpdate() { - int currentStep = m_GlobalTracker->GetCurrentStep(); - mitk::ProgressBar::GetInstance()->Progress(currentStep-m_LastStep); - UpdateTrackingStatus(); - GenerateFiberBundle(false); - m_LastStep = currentStep; + int currentStep = m_GlobalTracker->GetCurrentStep(); + mitk::ProgressBar::GetInstance()->Progress(currentStep-m_LastStep); + UpdateTrackingStatus(); + GenerateFiberBundle(false); + m_LastStep = currentStep; } // tell global tractography filter to stop after current step void QmitkGibbsTrackingView::StopGibbsTracking() { - if (m_GlobalTracker.IsNull()) - return; + if (m_GlobalTracker.IsNull()) + return; - //mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetSteps()-m_LastStep+1); - m_GlobalTracker->SetAbortTracking(true); - m_Controls->m_TrackingStop->setEnabled(false); - m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); + //mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetSteps()-m_LastStep+1); + m_GlobalTracker->SetAbortTracking(true); + m_Controls->m_TrackingStop->setEnabled(false); + m_Controls->m_TrackingStop->setText("Stopping Tractography ..."); } // update gui elements and generate fiber bundle after tracking is finished void QmitkGibbsTrackingView::AfterThread() { - m_ThreadIsRunning = false; - m_TrackingTimer->stop(); - - mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetSteps()-m_LastStep+1); - UpdateGUI(); - UpdateTrackingStatus(); - - if(m_Controls->m_ParticleWeightSlider->value()==0) - { - m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); - m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); - } - if(m_Controls->m_ParticleWidthSlider->value()==0) - { - m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); - m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); - } - if(m_Controls->m_ParticleLengthSlider->value()==0) - { - m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); - m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); - } - - GenerateFiberBundle(true); - m_FiberBundleNode = NULL; + m_ThreadIsRunning = false; + m_TrackingTimer->stop(); + + mitk::ProgressBar::GetInstance()->Progress(m_GlobalTracker->GetSteps()-m_LastStep+1); + UpdateGUI(); + UpdateTrackingStatus(); + + if(m_Controls->m_ParticleWeightSlider->value()==0) + { + m_Controls->m_ParticleWeightLabel->setText(QString::number(m_GlobalTracker->GetParticleWeight())); + m_Controls->m_ParticleWeightSlider->setValue(m_GlobalTracker->GetParticleWeight()*10000); + } + if(m_Controls->m_ParticleWidthSlider->value()==0) + { + m_Controls->m_ParticleWidthLabel->setText(QString::number(m_GlobalTracker->GetParticleWidth())); + m_Controls->m_ParticleWidthSlider->setValue(m_GlobalTracker->GetParticleWidth()*10); + } + if(m_Controls->m_ParticleLengthSlider->value()==0) + { + m_Controls->m_ParticleLengthLabel->setText(QString::number(m_GlobalTracker->GetParticleLength())); + m_Controls->m_ParticleLengthSlider->setValue(m_GlobalTracker->GetParticleLength()*10); + } + + GenerateFiberBundle(true); + m_FiberBundleNode = NULL; } // start tracking timer and update gui elements before tracking is started void QmitkGibbsTrackingView::BeforeThread() { - m_ThreadIsRunning = true; - m_TrackingTime = QTime::currentTime(); - m_ElapsedTime = 0; - m_TrackingTimer->start(1000); - m_LastStep = 0; + m_ThreadIsRunning = true; + m_TrackingTime = QTime::currentTime(); + m_ElapsedTime = 0; + m_TrackingTimer->start(1000); + m_LastStep = 0; - UpdateGUI(); + UpdateGUI(); } // setup gui elements and signal/slot connections void QmitkGibbsTrackingView::CreateQtPartControl( QWidget *parent ) { - // build up qt view, unless already done - if ( !m_Controls ) - { - // create GUI widgets from the Qt Designer's .ui file - m_Controls = new Ui::QmitkGibbsTrackingViewControls; - m_Controls->setupUi( parent ); - - AdvancedSettings(); - - connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); - connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); - connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); - connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); - connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); - connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); - connect( m_Controls->m_IterationsSlider, SIGNAL(valueChanged(int)), this, SLOT(SetIterations(int)) ); - connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); - connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); - connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); - connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); - connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); - connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); - connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); - connect( m_Controls->m_CurvatureThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(SetCurvatureThreshold(int)) ); - connect( m_Controls->m_OutputFileButton, SIGNAL(clicked()), this, SLOT(SetOutputFile()) ); - } + // build up qt view, unless already done + if ( !m_Controls ) + { + // create GUI widgets from the Qt Designer's .ui file + m_Controls = new Ui::QmitkGibbsTrackingViewControls; + m_Controls->setupUi( parent ); + + AdvancedSettings(); + + connect( m_TrackingTimer, SIGNAL(timeout()), this, SLOT(TimerUpdate()) ); + connect( m_Controls->m_TrackingStop, SIGNAL(clicked()), this, SLOT(StopGibbsTracking()) ); + connect( m_Controls->m_TrackingStart, SIGNAL(clicked()), this, SLOT(StartGibbsTracking()) ); + connect( m_Controls->m_AdvancedSettingsCheckbox, SIGNAL(clicked()), this, SLOT(AdvancedSettings()) ); + connect( m_Controls->m_SaveTrackingParameters, SIGNAL(clicked()), this, SLOT(SaveTrackingParameters()) ); + connect( m_Controls->m_LoadTrackingParameters, SIGNAL(clicked()), this, SLOT(LoadTrackingParameters()) ); + connect( m_Controls->m_IterationsSlider, SIGNAL(valueChanged(int)), this, SLOT(SetIterations(int)) ); + connect( m_Controls->m_ParticleWidthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWidth(int)) ); + connect( m_Controls->m_ParticleLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleLength(int)) ); + connect( m_Controls->m_InExBalanceSlider, SIGNAL(valueChanged(int)), this, SLOT(SetInExBalance(int)) ); + connect( m_Controls->m_FiberLengthSlider, SIGNAL(valueChanged(int)), this, SLOT(SetFiberLength(int)) ); + connect( m_Controls->m_ParticleWeightSlider, SIGNAL(valueChanged(int)), this, SLOT(SetParticleWeight(int)) ); + connect( m_Controls->m_StartTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetStartTemp(int)) ); + connect( m_Controls->m_EndTempSlider, SIGNAL(valueChanged(int)), this, SLOT(SetEndTemp(int)) ); + connect( m_Controls->m_CurvatureThresholdSlider, SIGNAL(valueChanged(int)), this, SLOT(SetCurvatureThreshold(int)) ); + connect( m_Controls->m_OutputFileButton, SIGNAL(clicked()), this, SLOT(SetOutputFile()) ); + } } void QmitkGibbsTrackingView::SetInExBalance(int value) { - m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); + m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); } void QmitkGibbsTrackingView::SetFiberLength(int value) { - m_Controls->m_FiberLengthLabel->setText(QString::number(value)+"mm"); + m_Controls->m_FiberLengthLabel->setText(QString::number(value)+"mm"); } void QmitkGibbsTrackingView::SetParticleWeight(int value) { - if (value>0) - m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); - else - m_Controls->m_ParticleWeightLabel->setText("auto"); + if (value>0) + m_Controls->m_ParticleWeightLabel->setText(QString::number((float)value/10000)); + else + m_Controls->m_ParticleWeightLabel->setText("auto"); } void QmitkGibbsTrackingView::SetStartTemp(int value) { - m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); + m_Controls->m_StartTempLabel->setText(QString::number((float)value/100)); } void QmitkGibbsTrackingView::SetEndTemp(int value) { - m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); + m_Controls->m_EndTempLabel->setText(QString::number((float)value/10000)); } void QmitkGibbsTrackingView::SetParticleWidth(int value) { - if (value>0) - m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); - else - m_Controls->m_ParticleWidthLabel->setText("auto"); + if (value>0) + m_Controls->m_ParticleWidthLabel->setText(QString::number((float)value/10)+" mm"); + else + m_Controls->m_ParticleWidthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetParticleLength(int value) { - if (value>0) - m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); - else - m_Controls->m_ParticleLengthLabel->setText("auto"); + if (value>0) + m_Controls->m_ParticleLengthLabel->setText(QString::number((float)value/10)+" mm"); + else + m_Controls->m_ParticleLengthLabel->setText("auto"); } void QmitkGibbsTrackingView::SetCurvatureThreshold(int value) { - m_Controls->m_CurvatureThresholdLabel->setText(QString::number(value)+"°"); + m_Controls->m_CurvatureThresholdLabel->setText(QString::number(value)+"°"); } void QmitkGibbsTrackingView::SetIterations(int value) { - switch(value) - { - case 0: - m_Controls->m_IterationsLabel->setText("Iterations: 1x10^4"); - m_Iterations = 10000; - break; - case 1: - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); - m_Iterations = 50000; - break; - case 2: - m_Controls->m_IterationsLabel->setText("Iterations: 1x10^5"); - m_Iterations = 100000; - break; - case 3: - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); - m_Iterations = 500000; - break; - case 4: - m_Controls->m_IterationsLabel->setText("Iterations: 1x10^6"); - m_Iterations = 1000000; - break; - case 5: - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); - m_Iterations = 5000000; - break; - case 6: - m_Controls->m_IterationsLabel->setText("Iterations: 1x10^7"); - m_Iterations = 10000000; - break; - case 7: - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); - m_Iterations = 50000000; - break; - case 8: - m_Controls->m_IterationsLabel->setText("Iterations: 1x10^8"); - m_Iterations = 100000000; - break; - case 9: - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); - m_Iterations = 500000000; - break; - case 10: - m_Controls->m_IterationsLabel->setText("Iterations: 1x10^9"); - m_Iterations = 1000000000; - break; - } + switch(value) + { + case 0: + m_Controls->m_IterationsLabel->setText("Iterations: 1x10^4"); + m_Iterations = 10000; + break; + case 1: + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); + m_Iterations = 50000; + break; + case 2: + m_Controls->m_IterationsLabel->setText("Iterations: 1x10^5"); + m_Iterations = 100000; + break; + case 3: + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); + m_Iterations = 500000; + break; + case 4: + m_Controls->m_IterationsLabel->setText("Iterations: 1x10^6"); + m_Iterations = 1000000; + break; + case 5: + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); + m_Iterations = 5000000; + break; + case 6: + m_Controls->m_IterationsLabel->setText("Iterations: 1x10^7"); + m_Iterations = 10000000; + break; + case 7: + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); + m_Iterations = 50000000; + break; + case 8: + m_Controls->m_IterationsLabel->setText("Iterations: 1x10^8"); + m_Iterations = 100000000; + break; + case 9: + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); + m_Iterations = 500000000; + break; + case 10: + m_Controls->m_IterationsLabel->setText("Iterations: 1x10^9"); + m_Iterations = 1000000000; + break; + } } void QmitkGibbsTrackingView::StdMultiWidgetAvailable(QmitkStdMultiWidget &stdMultiWidget) { - m_MultiWidget = &stdMultiWidget; + m_MultiWidget = &stdMultiWidget; } void QmitkGibbsTrackingView::StdMultiWidgetNotAvailable() { - m_MultiWidget = NULL; + m_MultiWidget = NULL; } // called if datamanager selection changes void QmitkGibbsTrackingView::OnSelectionChanged( std::vector nodes ) { - if (m_ThreadIsRunning) - return; - - m_QBallImageNode = NULL; - m_MaskImageNode = NULL; + if (m_ThreadIsRunning) + return; - // iterate all selected objects - for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) - { - mitk::DataNode::Pointer node = *it; + m_ImageNode = NULL; + m_MaskImageNode = NULL; - if( node.IsNotNull() && dynamic_cast(node->GetData()) ) - m_QBallImageNode = node; - else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + // iterate all selected objects + for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { - bool isBinary = false; - node->GetPropertyValue("binary", isBinary); - if (isBinary) - m_MaskImageNode = node; + mitk::DataNode::Pointer node = *it; + + if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + m_ImageNode = node; + else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + m_ImageNode = node; + else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) + { + bool isBinary = false; + node->GetPropertyValue("binary", isBinary); + if (isBinary) + m_MaskImageNode = node; + } } - } - UpdateGUI(); + UpdateGUI(); } // update gui elements displaying trackings status void QmitkGibbsTrackingView::UpdateTrackingStatus() { - if (m_GlobalTracker.IsNull()) - return; + if (m_GlobalTracker.IsNull()) + return; - m_ElapsedTime += m_TrackingTime.elapsed()/1000; - m_TrackingTime.restart(); - unsigned long hours = m_ElapsedTime/3600; - unsigned long minutes = (m_ElapsedTime%3600)/60; - unsigned long seconds = m_ElapsedTime%60; + m_ElapsedTime += m_TrackingTime.elapsed()/1000; + m_TrackingTime.restart(); + unsigned long hours = m_ElapsedTime/3600; + unsigned long minutes = (m_ElapsedTime%3600)/60; + unsigned long seconds = m_ElapsedTime%60; - m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); + m_Controls->m_ProposalAcceptance->setText(QString::number(m_GlobalTracker->GetProposalAcceptance()*100)+"%"); - m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); - m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); - m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); - m_Controls->m_CurrentStepLabel->setText( QString::number(100*(float)m_GlobalTracker->GetCurrentStep()/m_GlobalTracker->GetSteps())+"%" ); - m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); + m_Controls->m_TrackingTimeLabel->setText( QString::number(hours)+QString("h ")+QString::number(minutes)+QString("m ")+QString::number(seconds)+QString("s") ); + m_Controls->m_NumConnectionsLabel->setText( QString::number(m_GlobalTracker->GetNumConnections()) ); + m_Controls->m_NumParticlesLabel->setText( QString::number(m_GlobalTracker->GetNumParticles()) ); + m_Controls->m_CurrentStepLabel->setText( QString::number(100*(float)(m_GlobalTracker->GetCurrentStep()-1)/m_GlobalTracker->GetSteps())+"%" ); + m_Controls->m_AcceptedFibersLabel->setText( QString::number(m_GlobalTracker->GetNumAcceptedFibers()) ); } // update gui elements (enable/disable elements and set tooltips) void QmitkGibbsTrackingView::UpdateGUI() { - if (m_QBallImageNode.IsNotNull()) - m_Controls->m_QballImageLabel->setText(m_QBallImageNode->GetName().c_str()); - else - m_Controls->m_QballImageLabel->setText("-"); - if (m_MaskImageNode.IsNotNull()) - m_Controls->m_MaskImageLabel->setText(m_MaskImageNode->GetName().c_str()); - else - m_Controls->m_MaskImageLabel->setText("-"); - - if (!m_ThreadIsRunning && m_QBallImageNode.IsNotNull()) - { - m_Controls->m_TrackingStop->setEnabled(false); - m_Controls->m_TrackingStart->setEnabled(true); - m_Controls->m_LoadTrackingParameters->setEnabled(true); - m_Controls->m_IterationsSlider->setEnabled(true); - m_Controls->m_AdvancedFrame->setEnabled(true); - m_Controls->m_TrackingStop->setText("Stop Tractography"); - m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); - m_Controls->m_TrackingStop->setToolTip(""); - } - else if (!m_ThreadIsRunning) - { - m_Controls->m_TrackingStop->setEnabled(false); - m_Controls->m_TrackingStart->setEnabled(false); - m_Controls->m_LoadTrackingParameters->setEnabled(true); - m_Controls->m_IterationsSlider->setEnabled(true); - m_Controls->m_AdvancedFrame->setEnabled(true); - m_Controls->m_TrackingStop->setText("Stop Tractography"); - m_Controls->m_TrackingStart->setToolTip("No Q-Ball image selected."); - m_Controls->m_TrackingStop->setToolTip(""); - } - else - { - m_Controls->m_TrackingStop->setEnabled(true); - m_Controls->m_TrackingStart->setEnabled(false); - m_Controls->m_LoadTrackingParameters->setEnabled(false); - m_Controls->m_IterationsSlider->setEnabled(false); - m_Controls->m_AdvancedFrame->setEnabled(false); - m_Controls->m_AdvancedFrame->setVisible(false); - m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); - m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); - m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); - } + if (m_ImageNode.IsNotNull()) + m_Controls->m_QballImageLabel->setText(m_ImageNode->GetName().c_str()); + else + m_Controls->m_QballImageLabel->setText("-"); + if (m_MaskImageNode.IsNotNull()) + m_Controls->m_MaskImageLabel->setText(m_MaskImageNode->GetName().c_str()); + else + m_Controls->m_MaskImageLabel->setText("-"); + + if (!m_ThreadIsRunning && m_ImageNode.IsNotNull()) + { + m_Controls->m_TrackingStop->setEnabled(false); + m_Controls->m_TrackingStart->setEnabled(true); + m_Controls->m_LoadTrackingParameters->setEnabled(true); + m_Controls->m_IterationsSlider->setEnabled(true); + m_Controls->m_AdvancedFrame->setEnabled(true); + m_Controls->m_TrackingStop->setText("Stop Tractography"); + m_Controls->m_TrackingStart->setToolTip("Start tractography. No further change of parameters possible."); + m_Controls->m_TrackingStop->setToolTip(""); + } + else if (!m_ThreadIsRunning) + { + m_Controls->m_TrackingStop->setEnabled(false); + m_Controls->m_TrackingStart->setEnabled(false); + m_Controls->m_LoadTrackingParameters->setEnabled(true); + m_Controls->m_IterationsSlider->setEnabled(true); + m_Controls->m_AdvancedFrame->setEnabled(true); + m_Controls->m_TrackingStop->setText("Stop Tractography"); + m_Controls->m_TrackingStart->setToolTip("No Q-Ball image selected."); + m_Controls->m_TrackingStop->setToolTip(""); + } + else + { + m_Controls->m_TrackingStop->setEnabled(true); + m_Controls->m_TrackingStart->setEnabled(false); + m_Controls->m_LoadTrackingParameters->setEnabled(false); + m_Controls->m_IterationsSlider->setEnabled(false); + m_Controls->m_AdvancedFrame->setEnabled(false); + m_Controls->m_AdvancedFrame->setVisible(false); + m_Controls->m_AdvancedSettingsCheckbox->setChecked(false); + m_Controls->m_TrackingStart->setToolTip("Tracking in progress."); + m_Controls->m_TrackingStop->setToolTip("Stop tracking and display results."); + } } // show/hide advanced settings frame void QmitkGibbsTrackingView::AdvancedSettings() { - m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); + m_Controls->m_AdvancedFrame->setVisible(m_Controls->m_AdvancedSettingsCheckbox->isChecked()); } // set mask image data node void QmitkGibbsTrackingView::SetMask() { - std::vector nodes = GetDataManagerSelection(); - if (nodes.empty()) - { - m_MaskImageNode = NULL; - m_Controls->m_MaskImageLabel->setText("-"); - return; - } - - for( std::vector::iterator it = nodes.begin(); - it != nodes.end(); - ++it ) - { - mitk::DataNode::Pointer node = *it; + std::vector nodes = GetDataManagerSelection(); + if (nodes.empty()) + { + m_MaskImageNode = NULL; + m_Controls->m_MaskImageLabel->setText("-"); + return; + } - if (node.IsNotNull() && dynamic_cast(node->GetData())) + for( std::vector::iterator it = nodes.begin(); + it != nodes.end(); + ++it ) { - m_MaskImageNode = node; - m_Controls->m_MaskImageLabel->setText(node->GetName().c_str()); - return; + mitk::DataNode::Pointer node = *it; + + if (node.IsNotNull() && dynamic_cast(node->GetData())) + { + m_MaskImageNode = node; + m_Controls->m_MaskImageLabel->setText(node->GetName().c_str()); + return; + } } - } } // cast image to float template void QmitkGibbsTrackingView::CastToFloat(InputImageType* image, mitk::Image::Pointer outImage) { - typedef itk::CastImageFilter ItkCastFilter; - typename ItkCastFilter::Pointer itkCaster = ItkCastFilter::New(); - itkCaster->SetInput(image); - itkCaster->Update(); - outImage->InitializeByItk(itkCaster->GetOutput()); - outImage->SetVolume(itkCaster->GetOutput()->GetBufferPointer()); + typedef itk::CastImageFilter ItkCastFilter; + typename ItkCastFilter::Pointer itkCaster = ItkCastFilter::New(); + itkCaster->SetInput(image); + itkCaster->Update(); + outImage->InitializeByItk(itkCaster->GetOutput()); + outImage->SetVolume(itkCaster->GetOutput()->GetBufferPointer()); } // check for mask and qbi and start tracking thread void QmitkGibbsTrackingView::StartGibbsTracking() { - if(m_ThreadIsRunning) - { - MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; - return; - } - - if (m_QBallImageNode.IsNull()) - { - // Nothing selected. Inform the user and return - QMessageBox::information( NULL, "Warning", "Please load and select a qball image before starting image processing."); - return; - } - - // a node itself is not very useful, we need its data item (the image) - mitk::BaseData* data = m_QBallImageNode->GetData(); - if (!data) - return; - - // test if this data item is an image or not (could also be a surface or something totally different) - m_QBallImage = dynamic_cast( data ); - if (m_QBallImage.IsNull()) - return; - - // cast qbi to itk - m_ItkQBallImage = ItkQBallImgType::New(); - mitk::CastToItkImage(m_QBallImage, m_ItkQBallImage); - - // mask image found? - // catch exceptions thrown by the itkAccess macros - try{ - if(m_MaskImageNode.IsNotNull()) + if(m_ThreadIsRunning) { - m_MaskImage = 0; - if (dynamic_cast(m_MaskImageNode->GetData())) + MITK_WARN("QmitkGibbsTrackingView")<<"Thread already running!"; + return; + } - mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), - m_MaskImage); + if (m_ImageNode.IsNull()) + { + QMessageBox::information( NULL, "Warning", "Please load and select a qball image before starting image processing."); + return; } - } - catch(...) - { - QMessageBox::warning(NULL, "Warning", "Incompatible mask image chosen. Processing without masking."); - //reset mask image + + if (dynamic_cast(m_ImageNode->GetData())) + m_QBallImage = dynamic_cast(m_ImageNode->GetData()); + else if (dynamic_cast(m_ImageNode->GetData())) + m_TensorImage = dynamic_cast(m_ImageNode->GetData()); + + if (m_QBallImage.IsNull() && m_TensorImage.IsNull()) + return; + + // cast qbi to itk + m_ItkTensorImage = NULL; + m_ItkQBallImage = NULL; m_MaskImage = NULL; - } - unsigned int steps = m_Iterations/10000; - if (steps<10) - steps = 10; + if (m_QBallImage.IsNotNull()) + { + m_ItkQBallImage = ItkQBallImgType::New(); + mitk::CastToItkImage(m_QBallImage, m_ItkQBallImage); + } + else + { + m_ItkTensorImage = ItkTensorImage::New(); + mitk::CastToItkImage(m_TensorImage, m_ItkTensorImage); + } + + // mask image found? + // catch exceptions thrown by the itkAccess macros + try{ + if(m_MaskImageNode.IsNotNull()) + { + if (dynamic_cast(m_MaskImageNode->GetData())) + mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); + } + } + catch(...){}; + + unsigned int steps = m_Iterations/10000; + if (steps<10) + steps = 10; - m_LastStep = 1; - mitk::ProgressBar::GetInstance()->AddStepsToDo(steps); + m_LastStep = 1; + mitk::ProgressBar::GetInstance()->AddStepsToDo(steps); - // start worker thread - m_TrackingThread.start(QThread::LowestPriority); + // start worker thread + m_TrackingThread.start(QThread::LowestPriority); } // generate mitkFiberBundle from tracking filter output void QmitkGibbsTrackingView::GenerateFiberBundle(bool smoothFibers) { - if (m_GlobalTracker.IsNull() || (!(m_Controls->m_VisualizationCheckbox->isChecked() || m_Controls->m_VisualizeOnceButton->isChecked()) && m_ThreadIsRunning)) - return; - - if (m_Controls->m_VisualizeOnceButton->isChecked()) - m_Controls->m_VisualizeOnceButton->setChecked(false); - - vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); - if ( fiberBundle->GetNumberOfLines()==0 ) - return; - m_FiberBundle = mitk::FiberBundleX::New(fiberBundle); - -// if (smoothFibers) -// m_FiberBundle->DoFiberSmoothing(10); - - if (m_FiberBundleNode.IsNotNull()){ - GetDefaultDataStorage()->Remove(m_FiberBundleNode); - m_FiberBundleNode = 0; - } - m_FiberBundleNode = mitk::DataNode::New(); - m_FiberBundleNode->SetData(m_FiberBundle); - - QString name(m_QBallImageNode->GetName().c_str()); - name += "_FiberBundle"; - m_FiberBundleNode->SetName(name.toStdString()); - m_FiberBundleNode->SetVisibility(true); - - if (!m_OutputFileName.isEmpty()) - { - QString filename = m_OutputFileName; - mitk::FiberBundleXWriter::Pointer writer = mitk::FiberBundleXWriter::New(); - writer->SetFileName(filename.toStdString()); - writer->SetInputFiberBundleX(m_FiberBundle.GetPointer()); - try - { - writer->Update(); - QMessageBox::information(NULL, "Fiber bundle saved to", filename); + if (m_GlobalTracker.IsNull() || (!(m_Controls->m_VisualizationCheckbox->isChecked() || m_Controls->m_VisualizeOnceButton->isChecked()) && m_ThreadIsRunning)) + return; + + if (m_Controls->m_VisualizeOnceButton->isChecked()) + m_Controls->m_VisualizeOnceButton->setChecked(false); + + vtkSmartPointer fiberBundle = m_GlobalTracker->GetFiberBundle(); + if ( fiberBundle->GetNumberOfLines()==0 ) + return; + m_FiberBundle = mitk::FiberBundleX::New(fiberBundle); + + // if (smoothFibers) + // m_FiberBundle->DoFiberSmoothing(10); + + if (m_FiberBundleNode.IsNotNull()){ + GetDefaultDataStorage()->Remove(m_FiberBundleNode); + m_FiberBundleNode = 0; } - catch (itk::ExceptionObject &ex) - { - QMessageBox::information(NULL, "Fiber bundle could not be saved", QString("%1\n%2\n%3\n%4\n%5\n%6").arg(ex.GetNameOfClass()).arg(ex.GetFile()).arg(ex.GetLine()).arg(ex.GetLocation()).arg(ex.what()).arg(ex.GetDescription())); + m_FiberBundleNode = mitk::DataNode::New(); + m_FiberBundleNode->SetData(m_FiberBundle); + + QString name(m_ImageNode->GetName().c_str()); + name += "_FiberBundle"; + m_FiberBundleNode->SetName(name.toStdString()); + m_FiberBundleNode->SetVisibility(true); - if(m_QBallImageNode.IsNull()) - GetDataStorage()->Add(m_FiberBundleNode); - else - GetDataStorage()->Add(m_FiberBundleNode, m_QBallImageNode); + if (!m_OutputFileName.isEmpty()) + { + QString filename = m_OutputFileName; + mitk::FiberBundleXWriter::Pointer writer = mitk::FiberBundleXWriter::New(); + writer->SetFileName(filename.toStdString()); + writer->SetInputFiberBundleX(m_FiberBundle.GetPointer()); + try + { + writer->Update(); + QMessageBox::information(NULL, "Fiber bundle saved to", filename); + } + catch (itk::ExceptionObject &ex) + { + QMessageBox::information(NULL, "Fiber bundle could not be saved", QString("%1\n%2\n%3\n%4\n%5\n%6").arg(ex.GetNameOfClass()).arg(ex.GetFile()).arg(ex.GetLine()).arg(ex.GetLocation()).arg(ex.what()).arg(ex.GetDescription())); + + if(m_ImageNode.IsNull()) + GetDataStorage()->Add(m_FiberBundleNode); + else + GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); + } + } + else { + if(m_ImageNode.IsNull()) + GetDataStorage()->Add(m_FiberBundleNode); + else + GetDataStorage()->Add(m_FiberBundleNode, m_ImageNode); } - } - else { - if(m_QBallImageNode.IsNull()) - GetDataStorage()->Add(m_FiberBundleNode); - else - GetDataStorage()->Add(m_FiberBundleNode, m_QBallImageNode); - } } void QmitkGibbsTrackingView::SetOutputFile() { - // SELECT FOLDER DIALOG - m_OutputFileName = QFileDialog::getSaveFileName(0, - tr("Set file name"), - QDir::currentPath()+"/FiberBundle.fib", - tr("Fiber Bundle (*.fib)") ); - if (m_OutputFileName.isEmpty()) - m_Controls->m_OutputFileLabel->setText("N/A"); - else - m_Controls->m_OutputFileLabel->setText(m_OutputFileName); + // SELECT FOLDER DIALOG + m_OutputFileName = QFileDialog::getSaveFileName(0, + tr("Set file name"), + QDir::currentPath()+"/FiberBundle.fib", + tr("Fiber Bundle (*.fib)") ); + if (m_OutputFileName.isEmpty()) + m_Controls->m_OutputFileLabel->setText("N/A"); + else + m_Controls->m_OutputFileLabel->setText(m_OutputFileName); } // save current tracking paramters as xml file (.gtp) void QmitkGibbsTrackingView::SaveTrackingParameters() { - TiXmlDocument documentXML; - TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); - documentXML.LinkEndChild( declXML ); - - TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); - mainXML->SetAttribute("file_version", "0.1"); - documentXML.LinkEndChild(mainXML); - - TiXmlElement* paramXML = new TiXmlElement("parameter_set"); - paramXML->SetAttribute("iterations", QString::number(m_Iterations).toStdString()); - paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); - paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); - paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); - paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); - paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); - paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); - paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); - paramXML->SetAttribute("curvature_threshold", QString::number(m_Controls->m_CurvatureThresholdSlider->value()).toStdString()); - mainXML->LinkEndChild(paramXML); - QString filename = QFileDialog::getSaveFileName( - 0, - tr("Save Parameters"), - QDir::currentPath()+"/param.gtp", - tr("Global Tracking Parameters (*.gtp)") ); - - if(filename.isEmpty() || filename.isNull()) - return; - if(!filename.endsWith(".gtp")) - filename += ".gtp"; - documentXML.SaveFile( filename.toStdString() ); + TiXmlDocument documentXML; + TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); + documentXML.LinkEndChild( declXML ); + + TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); + mainXML->SetAttribute("file_version", "0.1"); + documentXML.LinkEndChild(mainXML); + + TiXmlElement* paramXML = new TiXmlElement("parameter_set"); + paramXML->SetAttribute("iterations", QString::number(m_Iterations).toStdString()); + paramXML->SetAttribute("particle_length", QString::number((float)m_Controls->m_ParticleLengthSlider->value()/10).toStdString()); + paramXML->SetAttribute("particle_width", QString::number((float)m_Controls->m_ParticleWidthSlider->value()/10).toStdString()); + paramXML->SetAttribute("particle_weight", QString::number((float)m_Controls->m_ParticleWeightSlider->value()/10000).toStdString()); + paramXML->SetAttribute("temp_start", QString::number((float)m_Controls->m_StartTempSlider->value()/100).toStdString()); + paramXML->SetAttribute("temp_end", QString::number((float)m_Controls->m_EndTempSlider->value()/10000).toStdString()); + paramXML->SetAttribute("inexbalance", QString::number((float)m_Controls->m_InExBalanceSlider->value()/10).toStdString()); + paramXML->SetAttribute("fiber_length", QString::number(m_Controls->m_FiberLengthSlider->value()).toStdString()); + paramXML->SetAttribute("curvature_threshold", QString::number(m_Controls->m_CurvatureThresholdSlider->value()).toStdString()); + mainXML->LinkEndChild(paramXML); + QString filename = QFileDialog::getSaveFileName( + 0, + tr("Save Parameters"), + QDir::currentPath()+"/param.gtp", + tr("Global Tracking Parameters (*.gtp)") ); + + if(filename.isEmpty() || filename.isNull()) + return; + if(!filename.endsWith(".gtp")) + filename += ".gtp"; + documentXML.SaveFile( filename.toStdString() ); } void QmitkGibbsTrackingView::UpdateIteraionsGUI(unsigned long iterations) { - switch(iterations) - { - case 10000: - m_Controls->m_IterationsSlider->setValue(0); - m_Controls->m_IterationsLabel->setText("Iterations: 10^4"); - break; - case 50000: - m_Controls->m_IterationsSlider->setValue(1); - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); - break; - case 100000: - m_Controls->m_IterationsSlider->setValue(2); - m_Controls->m_IterationsLabel->setText("Iterations: 10^5"); - break; - case 500000: - m_Controls->m_IterationsSlider->setValue(3); - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); - break; - case 1000000: - m_Controls->m_IterationsSlider->setValue(4); - m_Controls->m_IterationsLabel->setText("Iterations: 10^6"); - break; - case 5000000: - m_Controls->m_IterationsSlider->setValue(5); - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); - break; - case 10000000: - m_Controls->m_IterationsSlider->setValue(6); - m_Controls->m_IterationsLabel->setText("Iterations: 10^7"); - break; - case 50000000: - m_Controls->m_IterationsSlider->setValue(7); - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); - break; - case 100000000: - m_Controls->m_IterationsSlider->setValue(8); - m_Controls->m_IterationsLabel->setText("Iterations: 10^8"); - break; - case 500000000: - m_Controls->m_IterationsSlider->setValue(9); - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); - break; - case 1000000000: - m_Controls->m_IterationsSlider->setValue(10); - m_Controls->m_IterationsLabel->setText("Iterations: 10^9"); - break; - case 5000000000: - m_Controls->m_IterationsSlider->setValue(11); - m_Controls->m_IterationsLabel->setText("Iterations: 5x10^9"); - break; - } + switch(iterations) + { + case 10000: + m_Controls->m_IterationsSlider->setValue(0); + m_Controls->m_IterationsLabel->setText("Iterations: 10^4"); + break; + case 50000: + m_Controls->m_IterationsSlider->setValue(1); + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^4"); + break; + case 100000: + m_Controls->m_IterationsSlider->setValue(2); + m_Controls->m_IterationsLabel->setText("Iterations: 10^5"); + break; + case 500000: + m_Controls->m_IterationsSlider->setValue(3); + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^5"); + break; + case 1000000: + m_Controls->m_IterationsSlider->setValue(4); + m_Controls->m_IterationsLabel->setText("Iterations: 10^6"); + break; + case 5000000: + m_Controls->m_IterationsSlider->setValue(5); + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^6"); + break; + case 10000000: + m_Controls->m_IterationsSlider->setValue(6); + m_Controls->m_IterationsLabel->setText("Iterations: 10^7"); + break; + case 50000000: + m_Controls->m_IterationsSlider->setValue(7); + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^7"); + break; + case 100000000: + m_Controls->m_IterationsSlider->setValue(8); + m_Controls->m_IterationsLabel->setText("Iterations: 10^8"); + break; + case 500000000: + m_Controls->m_IterationsSlider->setValue(9); + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^8"); + break; + case 1000000000: + m_Controls->m_IterationsSlider->setValue(10); + m_Controls->m_IterationsLabel->setText("Iterations: 10^9"); + break; + case 5000000000: + m_Controls->m_IterationsSlider->setValue(11); + m_Controls->m_IterationsLabel->setText("Iterations: 5x10^9"); + break; + } } // load current tracking paramters from xml file (.gtp) void QmitkGibbsTrackingView::LoadTrackingParameters() { - QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); - if(filename.isEmpty() || filename.isNull()) - return; - - TiXmlDocument doc( filename.toStdString() ); - doc.LoadFile(); - - TiXmlHandle hDoc(&doc); - TiXmlElement* pElem; - TiXmlHandle hRoot(0); - - pElem = hDoc.FirstChildElement().Element(); - hRoot = TiXmlHandle(pElem); - pElem = hRoot.FirstChildElement("parameter_set").Element(); - - QString iterations(pElem->Attribute("iterations")); - m_Iterations = iterations.toULong(); - UpdateIteraionsGUI(m_Iterations); - - QString particleLength(pElem->Attribute("particle_length")); - float pLength = particleLength.toFloat(); - QString particleWidth(pElem->Attribute("particle_width")); - float pWidth = particleWidth.toFloat(); - - if (pLength==0) - m_Controls->m_ParticleLengthLabel->setText("auto"); - else - m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); - if (pWidth==0) - m_Controls->m_ParticleWidthLabel->setText("auto"); - else - m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); - - m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); - m_Controls->m_ParticleLengthSlider->setValue(pLength*10); - - QString partWeight(pElem->Attribute("particle_weight")); - m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); - m_Controls->m_ParticleWeightLabel->setText(partWeight); - - QString startTemp(pElem->Attribute("temp_start")); - m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); - m_Controls->m_StartTempLabel->setText(startTemp); - - QString endTemp(pElem->Attribute("temp_end")); - m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); - m_Controls->m_EndTempLabel->setText(endTemp); - - QString inExBalance(pElem->Attribute("inexbalance")); - m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); - m_Controls->m_InExBalanceLabel->setText(inExBalance); - - QString fiberLength(pElem->Attribute("fiber_length")); - m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); - m_Controls->m_FiberLengthLabel->setText(fiberLength+"mm"); - - QString curvThres(pElem->Attribute("curvature_threshold")); - m_Controls->m_CurvatureThresholdSlider->setValue(curvThres.toInt()); - m_Controls->m_CurvatureThresholdLabel->setText(curvThres+"°"); + QString filename = QFileDialog::getOpenFileName(0, tr("Load Parameters"), QDir::currentPath(), tr("Global Tracking Parameters (*.gtp)") ); + if(filename.isEmpty() || filename.isNull()) + return; + + TiXmlDocument doc( filename.toStdString() ); + doc.LoadFile(); + + TiXmlHandle hDoc(&doc); + TiXmlElement* pElem; + TiXmlHandle hRoot(0); + + pElem = hDoc.FirstChildElement().Element(); + hRoot = TiXmlHandle(pElem); + pElem = hRoot.FirstChildElement("parameter_set").Element(); + + QString iterations(pElem->Attribute("iterations")); + m_Iterations = iterations.toULong(); + UpdateIteraionsGUI(m_Iterations); + + QString particleLength(pElem->Attribute("particle_length")); + float pLength = particleLength.toFloat(); + QString particleWidth(pElem->Attribute("particle_width")); + float pWidth = particleWidth.toFloat(); + + if (pLength==0) + m_Controls->m_ParticleLengthLabel->setText("auto"); + else + m_Controls->m_ParticleLengthLabel->setText(particleLength+" mm"); + if (pWidth==0) + m_Controls->m_ParticleWidthLabel->setText("auto"); + else + m_Controls->m_ParticleWidthLabel->setText(particleWidth+" mm"); + + m_Controls->m_ParticleWidthSlider->setValue(pWidth*10); + m_Controls->m_ParticleLengthSlider->setValue(pLength*10); + + QString partWeight(pElem->Attribute("particle_weight")); + m_Controls->m_ParticleWeightSlider->setValue(partWeight.toFloat()*10000); + m_Controls->m_ParticleWeightLabel->setText(partWeight); + + QString startTemp(pElem->Attribute("temp_start")); + m_Controls->m_StartTempSlider->setValue(startTemp.toFloat()*100); + m_Controls->m_StartTempLabel->setText(startTemp); + + QString endTemp(pElem->Attribute("temp_end")); + m_Controls->m_EndTempSlider->setValue(endTemp.toFloat()*10000); + m_Controls->m_EndTempLabel->setText(endTemp); + + QString inExBalance(pElem->Attribute("inexbalance")); + m_Controls->m_InExBalanceSlider->setValue(inExBalance.toFloat()*10); + m_Controls->m_InExBalanceLabel->setText(inExBalance); + + QString fiberLength(pElem->Attribute("fiber_length")); + m_Controls->m_FiberLengthSlider->setValue(fiberLength.toInt()); + m_Controls->m_FiberLengthLabel->setText(fiberLength+"mm"); + + QString curvThres(pElem->Attribute("curvature_threshold")); + m_Controls->m_CurvatureThresholdSlider->setValue(curvThres.toInt()); + m_Controls->m_CurvatureThresholdLabel->setText(curvThres+"°"); } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h index d7019ced8b..4aca41d356 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.h @@ -1,165 +1,168 @@ /*=================================================================== 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 QmitkGibbsTrackingView_h #define QmitkGibbsTrackingView_h #include #include #include "ui_QmitkGibbsTrackingViewControls.h" #include #include #include #include #include #include #include +#include +#include class QmitkGibbsTrackingView; class QmitkTrackingWorker : public QObject { Q_OBJECT public: QmitkTrackingWorker(QmitkGibbsTrackingView* view); public slots: void run(); private: QmitkGibbsTrackingView* m_View; }; /*! \brief QmitkGibbsTrackingView \warning This application module is not yet documented. Use "svn blame/praise/annotate" and ask the author to provide basic documentation. \sa QmitkFunctionality \ingroup Functionalities */ typedef itk::Image< float, 3 > FloatImageType; namespace itk { template class GibbsTrackingFilter; } class QmitkGibbsTrackingView : public QmitkFunctionality { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: - typedef itk::Image MaskImgType; - - typedef itk::Vector OdfVectorType; - typedef itk::Image ItkQBallImgType; - - typedef itk::GibbsTrackingFilter< ItkQBallImgType > GibbsTrackingFilterType; + typedef itk::Image ItkFloatImageType; + typedef itk::Vector OdfVectorType; + typedef itk::Image ItkQBallImgType; + typedef itk::Image< itk::DiffusionTensor3D, 3 > ItkTensorImage; + typedef itk::GibbsTrackingFilter< ItkQBallImgType > GibbsTrackingFilterType; static const std::string VIEW_ID; QmitkGibbsTrackingView(); virtual ~QmitkGibbsTrackingView(); virtual void CreateQtPartControl(QWidget *parent); virtual void StdMultiWidgetAvailable (QmitkStdMultiWidget &stdMultiWidget); virtual void StdMultiWidgetNotAvailable(); signals: protected slots: void StartGibbsTracking(); void StopGibbsTracking(); void AfterThread(); void BeforeThread(); void TimerUpdate(); void SetMask(); void AdvancedSettings(); void SaveTrackingParameters(); void LoadTrackingParameters(); void SetIterations(int value); void SetParticleWidth(int value); void SetParticleLength(int value); void SetInExBalance(int value); void SetFiberLength(int value); void SetParticleWeight(int value); void SetStartTemp(int value); void SetEndTemp(int value); void SetCurvatureThreshold(int value); void SetOutputFile(); private: // Visualization & GUI void GenerateFiberBundle(bool smoothFibers); void UpdateGUI(); void UpdateTrackingStatus(); /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); template void CastToFloat(InputImageType* image, typename mitk::Image::Pointer outImage); void UpdateIteraionsGUI(unsigned long iterations); Ui::QmitkGibbsTrackingViewControls* m_Controls; QmitkStdMultiWidget* m_MultiWidget; // data objects mitk::FiberBundleX::Pointer m_FiberBundle; - MaskImgType::Pointer m_MaskImage; + ItkFloatImageType::Pointer m_MaskImage; + mitk::TensorImage::Pointer m_TensorImage; mitk::QBallImage::Pointer m_QBallImage; ItkQBallImgType::Pointer m_ItkQBallImage; + ItkTensorImage::Pointer m_ItkTensorImage; // data nodes - mitk::DataNode::Pointer m_QBallImageNode; + mitk::DataNode::Pointer m_ImageNode; mitk::DataNode::Pointer m_MaskImageNode; mitk::DataNode::Pointer m_FiberBundleNode; // flags etc. bool m_ThreadIsRunning; QTimer* m_TrackingTimer; QTime m_TrackingTime; unsigned long m_ElapsedTime; unsigned long m_Iterations; int m_LastStep; QString m_OutputFileName; // global tracker and friends itk::SmartPointer m_GlobalTracker; QmitkTrackingWorker m_TrackingWorker; QThread m_TrackingThread; friend class QmitkTrackingWorker; }; #endif // _QMITKGibbsTrackingVIEW_H_INCLUDED diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui index 2558f22cca..44fbab145e 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui @@ -1,1013 +1,1013 @@ QmitkGibbsTrackingViewControls 0 0 463 1011 0 0 0 0 QmitkTemplate 0 9 3 9 3 Data - Q-Ball Image: + Q-Ball/Tensor Image: Mandatory input - Mask Image: Optional input to limit the algorithms search space. - Parameters 0 Iterations: 10^7 Specify number of iterations for the tracking algorithm. 10 6 Qt::Horizontal QSlider::TicksBelow true Activate continuous visualization of intermediate results. Visualize Tractography true Visualize intermediate result. :/QmitkDiffusionImaging/Refresh_48.png:/QmitkDiffusionImaging/Refresh_48.png true Advanced Settings Output File: QFrame::NoFrame QFrame::Plain 0 0 0 Select output file name and folder. ... N/A true true QFrame::StyledPanel QFrame::Raised 9 0 9 0 4 auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 0 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter Particle Width: Only fibers longer than specified are accepted. 500 1 40 Qt::Horizontal QSlider::NoTicks Curvature Threshold: 0.001 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 40mm Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 0.1 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter Balance In/Ex Energy: IE Bias < 0 < EE Bias -50 50 1 Qt::Horizontal QSlider::NoTicks automatic estimation from gfa map and q-ball data. 0 1000 1 0 Qt::Horizontal true QSlider::NoTicks Min. Fiber Length: Particle Length: 1 99 1 10 Qt::Horizontal QSlider::NoTicks Particle Weight: Start Temperature: Allow only fiber curvature values smaller than the selected threshold. 180 1 45 Qt::Horizontal QSlider::NoTicks auto = 1.5 * min. spacing; l 100 1 Qt::Horizontal QSlider::NoTicks auto = 0.5 * min. spacing; sigma 100 1 Qt::Horizontal QSlider::NoTicks 45° Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter End Temperature: 1 100 1 10 Qt::Horizontal false false QSlider::NoTicks QFrame::NoFrame QFrame::Plain 0 0 0 true Save current parameters as xml (.gtp) Qt::LeftToRight Save Parameters :/qmitk/btnMoveDown.png:/qmitk/btnMoveDown.png true Load parameters from xml file (.gtp) Qt::LeftToRight Load Parameters :/qmitk/btnMoveUp.png:/qmitk/btnMoveUp.png false No Q-Ball image selected. Qt::LeftToRight Start Tractography :/qmitk/play.xpm:/qmitk/play.xpm false Qt::LeftToRight Stop Tractography :/qmitk/stop.xpm:/qmitk/stop.xpm Monitor Progress: - Will only be updated if tracking is visualized Will only be updated if tracking is visualized Accepted Fibers: Connections: Particles: Proposal Acceptance Rate: Tracking Time: Will only be updated if tracking is visualized - - - - - Qt::Vertical QSizePolicy::Expanding 0 0