diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp index 6412ef569a..27d0cf8f76 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp @@ -1,423 +1,443 @@ /*=================================================================== 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" using namespace mitk; EnergyComputer::EnergyComputer(MTRand* rgen, ItkQBallImgType* qballImage, SphereInterpolator *sp, ParticleGrid *particleGrid, float *mask, vnl_matrix_fixed rotMatrix) { + m_UseTrilinearInterpolation = true; m_ParticleGrid = particleGrid; m_RandGen = rgen; m_RotationMatrix = rotMatrix; m_ImageData = qballImage; m_SphereInterpolator = sp; m_MaskImageData = 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]; + 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) 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_NumActiveVoxels = 0; cumulspatprob[0] = 0; - for (int x = 1; x < m_Size[0];x++) - for (int y = 1; y < m_Size[1];y++) - for (int z = 1; z < m_Size[2];z++) + 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) { cumulspatprob[m_NumActiveVoxels+1] = cumulspatprob[m_NumActiveVoxels] + m_MaskImageData[idx]; spatidx[m_NumActiveVoxels] = idx; m_NumActiveVoxels++; } } for (int k = 0; k < m_NumActiveVoxels; k++) cumulspatprob[k] /= cumulspatprob[m_NumActiveVoxels]; fprintf(stderr,"EnergyComputer: %i active voxels found\n",m_NumActiveVoxels); } void EnergyComputer::setParameters(float pwei,float pwid,float chempot_connection, float length,float curv_hardthres, float inex_balance, float chempot2, float meanv) { this->chempot2 = chempot2; meanval_sq = meanv; eigencon_energy = chempot_connection; eigen_energy = 0; particle_weight = pwei; 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) particle_length_sq = length*length; curv_hard = curv_hardthres; float sigma_s = pwid; gamma_s = 1/(sigma_s*sigma_s); gamma_reg_s =1/(length*length/4); } void EnergyComputer::drawSpatPosition(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]) { rh = j; continue; } if (r > cumulspatprob[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()); } float EnergyComputer::SpatProb(vnl_vector_fixed R) { - int rx = int(R[0]/m_Spacing[0]); - int ry = int(R[1]/m_Spacing[1]); - int rz = int(R[2]/m_Spacing[2]); + 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)]; } else return 0; } float EnergyComputer::evaluateODF(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len) { 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++) { Rs = R + (N * len) * ((float)i/CU); - 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; - - xint = int(floor(Rx)); - yint = int(floor(Ry)); - zint = int(floor(Rz)); - - if (xint >= 0 && xint < m_Size[0]-1 && yint >= 0 && yint < m_Size[1]-1 && zint >= 0 && zint < m_Size[2]-1) + if (!m_UseTrilinearInterpolation) { - 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; - - 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; - - 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; - - 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; - - 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; - - 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; - - 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; - - 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; + 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]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]); + } + } + else + { + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + } } } Dn *= 1.0/(2*CU+1); return Dn; } 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 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); 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 EnergyComputer::computeInternalEnergy(Particle *dp) { float energy = eigen_energy; if (dp->pID != -1) energy += computeInternalEnergyConnection(dp,+1); if (dp->mID != -1) energy += computeInternalEnergyConnection(dp,-1); return energy; } float EnergyComputer::computeInternalEnergyConnection(Particle *p1,int ep1) { Particle *p2 = 0; int ep2; if (ep1 == 1) p2 = m_ParticleGrid->GetParticle(p1->pID); 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); } float EnergyComputer::computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) { if ((dot_product(p1->N,p2->N))*ep1*ep2 > -curv_hard) return -INFINITY; vnl_vector_fixed R1 = p1->R + (p1->N * (p1->len * ep1)); vnl_vector_fixed R2 = p2->R + (p2->N * (p2->len * ep2)); if ((R1-R2).squared_magnitude() > particle_length_sq) 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; 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 fc143f8974..6b73d389c9 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h @@ -1,104 +1,106 @@ /*=================================================================== 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; float chempot2; float meanval_sq; float gamma_s; float gamma_reg_s; float particle_weight; float ex_strength; float in_strength; float particle_length_sq; float curv_hard; EnergyComputer(MTRand* rgen, ItkQBallImgType* qballImage, SphereInterpolator *sp, ParticleGrid *pcon, float *mask, vnl_matrix_fixed rotMatrix); void setParameters(float pwei,float pwid,float chempot_connection, float length,float curv_hardthres, float inex_balance, float chempot2, float meanv); 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); float mbesseli0(float x); float mexp(float x); }; #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp index 2d639cbc6f..0446ee26b5 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp @@ -1,141 +1,141 @@ /*=================================================================== 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 "mitkFiberBuilder.h" using namespace mitk; FiberBuilder::FiberBuilder(ParticleGrid* grid, ItkFloatImageType* image) { m_Grid = grid; m_Image = image; m_FiberLength = 0; } FiberBuilder::~FiberBuilder() { } vtkSmartPointer FiberBuilder::iterate(int minFiberLength) { m_VtkCellArray = vtkSmartPointer::New(); m_VtkPoints = vtkSmartPointer::New(); int cur_label = 1; int numFibers = 0; m_FiberLength = 0; for (int k = 0; k < m_Grid->m_NumParticles;k++) { Particle *dp = m_Grid->GetParticle(k); if (dp->label == 0) { vtkSmartPointer container = vtkSmartPointer::New(); dp->label = cur_label; dp->numerator = 0; labelPredecessors(dp, container); labelSuccessors(dp, container); cur_label++; if(m_FiberLength >= minFiberLength) { m_VtkCellArray->InsertNextCell(container); numFibers++; } m_FiberLength = 0; } } for (int k = 0; k < m_Grid->m_NumParticles;k++) { Particle *dp = m_Grid->GetParticle(k); dp->inserted = false; dp->label = 0; } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(m_VtkPoints); fiberPolyData->SetLines(m_VtkCellArray); return fiberPolyData; } void FiberBuilder::AddPoint(Particle *dp, vtkSmartPointer container) { if (dp->inserted) return; dp->inserted = true; itk::ContinuousIndex index; - index[0] = dp->R[0]/m_Image->GetSpacing()[0]; - index[1] = dp->R[1]/m_Image->GetSpacing()[1]; - index[2] = dp->R[2]/m_Image->GetSpacing()[2]; + index[0] = dp->R[0]/m_Image->GetSpacing()[0]-0.5; + index[1] = dp->R[1]/m_Image->GetSpacing()[1]-0.5; + index[2] = dp->R[2]/m_Image->GetSpacing()[2]-0.5; itk::Point point; m_Image->TransformContinuousIndexToPhysicalPoint( index, point ); vtkIdType id = m_VtkPoints->InsertNextPoint(point.GetDataPointer()); container->GetPointIds()->InsertNextId(id); if(container->GetNumberOfPoints()>1) m_FiberLength += m_LastPoint.EuclideanDistanceTo(point); m_LastPoint = point; } void FiberBuilder::labelPredecessors(Particle *dp, vtkSmartPointer container) { if (dp->mID != -1 && dp->mID!=dp->ID) { if (dp->ID!=m_Grid->GetParticle(dp->mID)->pID) { if (dp->ID==m_Grid->GetParticle(dp->mID)->mID) { int tmp = m_Grid->GetParticle(dp->mID)->pID; m_Grid->GetParticle(dp->mID)->pID = m_Grid->GetParticle(dp->mID)->mID; m_Grid->GetParticle(dp->mID)->mID = tmp; } } if (m_Grid->GetParticle(dp->mID)->label == 0) { m_Grid->GetParticle(dp->mID)->label = dp->label; m_Grid->GetParticle(dp->mID)->numerator = dp->numerator-1; labelPredecessors(m_Grid->GetParticle(dp->mID), container); } } AddPoint(dp, container); } void FiberBuilder::labelSuccessors(Particle *dp, vtkSmartPointer container) { AddPoint(dp, container); if (dp->pID != -1 && dp->pID!=dp->ID) { if (dp->ID!=m_Grid->GetParticle(dp->pID)->mID) { if (dp->ID==m_Grid->GetParticle(dp->pID)->pID) { int tmp = m_Grid->GetParticle(dp->pID)->pID; m_Grid->GetParticle(dp->pID)->pID = m_Grid->GetParticle(dp->pID)->mID; m_Grid->GetParticle(dp->pID)->mID = tmp; } } if (m_Grid->GetParticle(dp->pID)->label == 0) { m_Grid->GetParticle(dp->pID)->label = dp->label; m_Grid->GetParticle(dp->pID)->numerator = dp->numerator+1; labelSuccessors(m_Grid->GetParticle(dp->pID), container); } } } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index be493149bb..e1d0916fea 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,397 +1,389 @@ /*=================================================================== 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 // MISC #include #include #include namespace itk{ template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::GibbsTrackingFilter(): m_TempStart(0.1), m_TempEnd(0.001), m_NumIt(500000), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), m_ChempotConnection(10), m_InexBalance(0), m_Chempot2(0.2), m_FiberLength(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_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() { 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)}; - // make sure image has enough slices - if (imgSize[1]<3 || imgSize[2]<3 || imgSize[3]<3) - { - MITK_INFO << "itkGibbsTrackingFilter: image size < 3 not supported"; - m_AbortTracking = true; - } - // 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; 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; } // perform mean subtraction on odfs typedef ImageRegionIterator< ItkQBallImageType > InputIteratorType; InputIteratorType it(qballImage, 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.IsNotNull()) - { - 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); - } - else + if(m_MaskImage.IsNull()) { - mask = 0; - maskImageSize[0] = imgSize[1]; - maskImageSize[1] = imgSize[2]; - maskImageSize[2] = imgSize[3]; + 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); } - int mask_oversamp_mult = maskImageSize[0]/imgSize[1]; + 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); // get paramters float minSpacing; if(spacing[0]m_NumIt) { MITK_INFO << "itkGibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } if (m_CurvatureHardThreshold < mitk::eps) m_CurvatureHardThreshold = 0; unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); MTRand randGen(1); srand(1); SphereInterpolator* interpolator = LoadSphereInterpolator(); MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; 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); MetropolisHastingsSampler* sampler = new MetropolisHastingsSampler(particleGrid, &encomp, &randGen, m_CurvatureHardThreshold); // main loop m_NumAcceptedFibers = 0; for( int 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); for (unsigned long i=0; iMakeProposal(); if (m_BuildFibers) { m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/m_NumIt; m_NumParticles = particleGrid->m_NumParticles; m_NumConnections = particleGrid->m_NumConnections; FiberBuilder fiberBuilder(particleGrid, m_MaskImage); m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); m_BuildFibers = false; } } } 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 > 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/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp index a0590f3e1b..656c2e207d 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp @@ -1,769 +1,750 @@ /*=================================================================== 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) { } 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(); } const std::string QmitkGibbsTrackingView::VIEW_ID = "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); } QmitkGibbsTrackingView::~QmitkGibbsTrackingView() { 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; } // tell global tractography filter to stop after current step void QmitkGibbsTrackingView::StopGibbsTracking() { 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 ..."); } // 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; } // 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; 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()) ); } } void QmitkGibbsTrackingView::SetInExBalance(int value) { m_Controls->m_InExBalanceLabel->setText(QString::number((float)value/10)); } void QmitkGibbsTrackingView::SetFiberLength(int value) { 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"); } void QmitkGibbsTrackingView::SetStartTemp(int value) { 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)); } 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"); } 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"); } void QmitkGibbsTrackingView::SetCurvatureThreshold(int 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; } } void QmitkGibbsTrackingView::StdMultiWidgetAvailable(QmitkStdMultiWidget &stdMultiWidget) { m_MultiWidget = &stdMultiWidget; } void QmitkGibbsTrackingView::StdMultiWidgetNotAvailable() { m_MultiWidget = NULL; } // called if datamanager selection changes void QmitkGibbsTrackingView::OnSelectionChanged( std::vector nodes ) { if (m_ThreadIsRunning) return; m_QBallImageNode = NULL; m_MaskImageNode = NULL; // iterate all selected objects for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast(node->GetData()) ) m_QBallImageNode = node; else if( node.IsNotNull() && dynamic_cast(node->GetData()) ) { bool isBinary = false; node->GetPropertyValue("binary", isBinary); if (isBinary) m_MaskImageNode = node; } } UpdateGUI(); } // update gui elements displaying trackings status void QmitkGibbsTrackingView::UpdateTrackingStatus() { 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_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()) ); } // 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."); } } // show/hide advanced settings frame void QmitkGibbsTrackingView::AdvancedSettings() { 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; 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()); } // 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()) { m_MaskImage = 0; if (dynamic_cast(m_MaskImageNode->GetData())) mitk::CastToItkImage(dynamic_cast(m_MaskImageNode->GetData()), m_MaskImage); } } catch(...) { QMessageBox::warning(NULL, "Warning", "Incompatible mask image chosen. Processing without masking."); //reset mask image m_MaskImage = NULL; } - - // if no mask image is selected generate it - if( m_MaskImage.IsNull() ) - { - m_MaskImage = MaskImgType::New(); - m_MaskImage->SetSpacing( m_ItkQBallImage->GetSpacing() ); // Set the image spacing - m_MaskImage->SetOrigin( m_ItkQBallImage->GetOrigin() ); // Set the image origin - m_MaskImage->SetDirection( m_ItkQBallImage->GetDirection() ); // Set the image direction - m_MaskImage->SetLargestPossibleRegion( m_ItkQBallImage->GetLargestPossibleRegion()); - m_MaskImage->SetBufferedRegion( m_ItkQBallImage->GetLargestPossibleRegion() ); - m_MaskImage->Allocate(); - - itk::ImageRegionIterator it (m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); - for (it = it.Begin(); !it.IsAtEnd(); ++it) - { - it.Set(1); - } - } - unsigned int steps = m_Iterations/10000; if (steps<10) steps = 10; m_LastStep = 1; mitk::ProgressBar::GetInstance()->AddStepsToDo(steps); // 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); } 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_QBallImageNode.IsNull()) GetDataStorage()->Add(m_FiberBundleNode); else GetDataStorage()->Add(m_FiberBundleNode, m_QBallImageNode); } } 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); } // 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() ); } 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; } } // 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+"°"); }