diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp index ec700a555a..d967001c01 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp @@ -1,339 +1,452 @@ /*=================================================================== 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, float *data, const int *dsz, double *cellsize, SphereInterpolator *sp, ParticleGrid *pcon, float *spimg, int spmult, vnl_matrix_fixed rotMatrix) +EnergyComputer::EnergyComputer(MTRand* rgen, ItkQBallImgType* data, const int *dsz, double *cellsize, SphereInterpolator *sp, ParticleGrid *pcon, float *spimg, int spmult, vnl_matrix_fixed rotMatrix) { mtrand = rgen; m_RotationMatrix = rotMatrix; m_QBallImageData = data; m_QBallImageSize = dsz; m_SphereInterpolator = sp; m_MaskImageData = spimg; nip = m_QBallImageSize[0]; w = m_QBallImageSize[1]; h = m_QBallImageSize[2]; d = m_QBallImageSize[3]; voxsize_w = cellsize[0]; voxsize_h = cellsize[1]; voxsize_d = cellsize[2]; w_sp = m_QBallImageSize[1]*spmult; h_sp = m_QBallImageSize[2]*spmult; d_sp = m_QBallImageSize[3]*spmult; voxsize_sp_w = cellsize[0]/spmult; voxsize_sp_h = cellsize[1]/spmult; voxsize_sp_d = cellsize[2]/spmult; fprintf(stderr,"Data size (voxels) : %i x %i x %i\n",w,h,d); fprintf(stderr,"voxel size: %f x %f x %f\n",voxsize_w,voxsize_h,voxsize_d); fprintf(stderr,"mask_oversamp_mult: %i\n",spmult); if (nip != sp->nverts) { fprintf(stderr,"EnergyComputer: error during init: data does not match with interpolation scheme\n"); } m_ParticleGrid = pcon; int totsz = w_sp*h_sp*d_sp; cumulspatprob = (float*) malloc(sizeof(float) * totsz); spatidx = (int*) malloc(sizeof(int) * totsz); if (cumulspatprob == 0 || spatidx == 0) { fprintf(stderr,"EnergyCOmputerBase: out of memory!\n"); return ; } scnt = 0; cumulspatprob[0] = 0; for (int x = 1; x < w_sp;x++) for (int y = 1; y < h_sp;y++) for (int z = 1; z < d_sp;z++) { int idx = x+(y+z*h_sp)*w_sp; if (m_MaskImageData[idx] > 0.5) { cumulspatprob[scnt+1] = cumulspatprob[scnt] + m_MaskImageData[idx]; spatidx[scnt] = idx; scnt++; } } for (int k = 0; k < scnt; k++) { cumulspatprob[k] /= cumulspatprob[scnt]; } fprintf(stderr,"#active voxels: %i (in mask units) \n",scnt); } 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 = mtrand->frand(); int j; int rl = 1; int rh = scnt; 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] = voxsize_sp_w*((float)(spatidx[rh-1] % w_sp) + mtrand->frand()); R[1] = voxsize_sp_h*((float)((spatidx[rh-1]/w_sp) % h_sp) + mtrand->frand()); R[2] = voxsize_sp_d*((float)(spatidx[rh-1]/(w_sp*h_sp)) + mtrand->frand()); } float EnergyComputer::SpatProb(vnl_vector_fixed R) { int rx = int(R[0]/voxsize_sp_w); int ry = int(R[1]/voxsize_sp_h); int rz = int(R[2]/voxsize_sp_d); if (rx >= 0 && rx < w_sp && ry >= 0 && ry < h_sp && rz >= 0 && rz < d_sp){ return m_MaskImageData[rx + w_sp* (ry + h_sp*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,spatindex; 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]/voxsize_w-0.5; float Ry = Rs[1]/voxsize_h-0.5; float Rz = Rs[2]/voxsize_d-0.5; xint = int(floor(Rx)); yint = int(floor(Ry)); zint = int(floor(Rz)); if (xint >= 0 && xint < w-1 && yint >= 0 && yint < h-1 && zint >= 0 && zint < d-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); - spatindex = (xint + w*(yint+h*zint)) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint; index[1] = yint; index[2] = zint; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(1-zfrac); - spatindex = (xint+1 + w*(yint+h*zint)) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint+1; index[1] = yint; index[2] = zint; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(1-zfrac); - spatindex = (xint + w*(yint+1+h*zint)) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint; index[1] = yint+1; index[2] = zint; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(1-yfrac)*(zfrac); - spatindex = (xint + w*(yint+h*(zint+1))) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint; index[1] = yint; index[2] = zint+1; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(yfrac)*(1-zfrac); - spatindex = (xint+1 + w*(yint+1+h*zint)) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint+1; index[1] = yint+1; index[2] = zint; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(zfrac); - spatindex = (xint + w*(yint+1+h*(zint+1))) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint; index[1] = yint+1; index[2] = zint+1; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(zfrac); - spatindex = (xint+1 + w*(yint+h*(zint+1))) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint+1; index[1] = yint; index[2] = zint+1; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(yfrac)*(zfrac); - spatindex = (xint+1 + w*(yint+1+h*(zint+1))) *nip; - Dn += (m_QBallImageData[spatindex + m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_QBallImageData[spatindex + m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; + index[0] = xint+1; index[1] = yint+1; index[2] = zint+1; + Dn += (m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + + m_QBallImageData->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + + m_QBallImageData->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 &cap, 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)*p->cap ; w = mexp(dpos*gamma_reg_s); Pn += w*bw; } } float energy = 0; energy += (2*(Dn/particle_weight-Sn) - (mbesseli0(1.0)+chempot2)*cap)*cap; 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->m_AddressContainer[p1->pID]; else p2 = m_ParticleGrid->m_AddressContainer[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 6a90719beb..4ccb73684a 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h @@ -1,97 +1,102 @@ /*=================================================================== 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 "mitkParticleGrid.h" #include "mitkSphereInterpolator.h" #include "MersenneTwister.h" #include #include using namespace mitk; class MitkDiffusionImaging_EXPORT EnergyComputer { public: + + typedef itk::Vector OdfVectorType; + typedef itk::Image ItkQBallImgType; + float eigen_energy; vnl_matrix_fixed m_RotationMatrix; - float *m_QBallImageData; + ItkQBallImgType* m_QBallImageData; const int *m_QBallImageSize; SphereInterpolator *m_SphereInterpolator; ParticleGrid *m_ParticleGrid; 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* mtrand; int nip; // number of data vertices on sphere float *m_MaskImageData; float *cumulspatprob; int *spatidx; int scnt; 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, float *data, const int *dsz, double *cellsize, SphereInterpolator *sp, ParticleGrid *pcon, float *spimg, int spmult, vnl_matrix_fixed rotMatrix); + EnergyComputer(MTRand* rgen, ItkQBallImgType* data, const int *dsz, double *cellsize, SphereInterpolator *sp, ParticleGrid *pcon, float *spimg, int spmult, 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 &cap, 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 095fa63b58..60c0bf17c3 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp @@ -1,147 +1,142 @@ /*=================================================================== 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(float *points, int numPoints, double spacing[], ItkQBallImgType::Pointer image) +FiberBuilder::FiberBuilder(ParticleGrid* grid, ItkFloatImageType* image) { + m_Grid = grid; + particles = m_Grid->m_Particles; + m_Image = image; m_FiberLength = 0; - m_ItkQBallImage = image; - particles = (Particle*) malloc(sizeof(Particle)*numPoints); - pcnt = numPoints; - attrcnt = 10; - for (int k = 0; k < numPoints; k++) - { - Particle *p = &(particles[k]); - p->R = vnl_vector_fixed(points[attrcnt*k]/spacing[0]-0.5, points[attrcnt*k+1]/spacing[1]-0.5,points[attrcnt*k+2]/spacing[2]-0.5); - p->N = vnl_vector_fixed(points[attrcnt*k+3],points[attrcnt*k+4],points[attrcnt*k+5]); - p->cap = points[attrcnt*k+6]; - p->len = points[attrcnt*k+7]; - p->mID = (int) points[attrcnt*k+8]; - p->pID = (int) points[attrcnt*k+9]; - p->ID = k; - p->label = 0; - } - m_VtkCellArray = vtkSmartPointer::New(); - m_VtkPoints = vtkSmartPointer::New(); } FiberBuilder::~FiberBuilder() { - free(particles); + } 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 < pcnt;k++) + for (int k = 0; k < m_Grid->m_NumParticles;k++) { Particle *dp = &(particles[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 = &(particles[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]; - index[1] = dp->R[1]; - index[2] = dp->R[2]; + 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]; itk::Point point; - m_ItkQBallImage->TransformContinuousIndexToPhysicalPoint( index, 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!=particles[dp->mID].pID) + if (dp->ID!=particles[m_Grid->Id2Index(dp->mID)].pID) { - if (dp->ID==particles[dp->mID].mID) + if (dp->ID==particles[m_Grid->Id2Index(dp->mID)].mID) { - int tmp = particles[dp->mID].pID; - particles[dp->mID].pID = particles[dp->mID].mID; - particles[dp->mID].mID = tmp; + int tmp = particles[m_Grid->Id2Index(dp->mID)].pID; + particles[m_Grid->Id2Index(dp->mID)].pID = particles[m_Grid->Id2Index(dp->mID)].mID; + particles[m_Grid->Id2Index(dp->mID)].mID = tmp; } } - if (particles[dp->mID].label == 0) + if (particles[m_Grid->Id2Index(dp->mID)].label == 0) { - particles[dp->mID].label = dp->label; - particles[dp->mID].numerator = dp->numerator-1; - labelPredecessors(&(particles[dp->mID]), container); + particles[m_Grid->Id2Index(dp->mID)].label = dp->label; + particles[m_Grid->Id2Index(dp->mID)].numerator = dp->numerator-1; + labelPredecessors(&(particles[m_Grid->Id2Index(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!=particles[dp->pID].mID) + if (dp->ID!=particles[m_Grid->Id2Index(dp->pID)].mID) { - if (dp->ID==particles[dp->pID].pID) + if (dp->ID==particles[m_Grid->Id2Index(dp->pID)].pID) { - int tmp = particles[dp->pID].pID; - particles[dp->pID].pID = particles[dp->pID].mID; - particles[dp->pID].mID = tmp; + int tmp = particles[m_Grid->Id2Index(dp->pID)].pID; + particles[m_Grid->Id2Index(dp->pID)].pID = particles[m_Grid->Id2Index(dp->pID)].mID; + particles[m_Grid->Id2Index(dp->pID)].mID = tmp; } } - if (particles[dp->pID].label == 0) + if (particles[m_Grid->Id2Index(dp->pID)].label == 0) { - particles[dp->pID].label = dp->label; - particles[dp->pID].numerator = dp->numerator+1; - labelSuccessors(&(particles[dp->pID]), container); + particles[m_Grid->Id2Index(dp->pID)].label = dp->label; + particles[m_Grid->Id2Index(dp->pID)].numerator = dp->numerator+1; + labelSuccessors(&(particles[m_Grid->Id2Index(dp->pID)]), container); } } } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h index 0f0593095a..a11485bf25 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h @@ -1,73 +1,70 @@ /*=================================================================== 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 _BUILDFIBRES #define _BUILDFIBRES #include #include #include "mitkParticleGrid.h" #include #include #include #include #include #include #include #include using namespace std; namespace mitk { class MitkDiffusionImaging_EXPORT FiberBuilder { public: - Particle *particles; - int pcnt; - int attrcnt; - typedef vector< Particle* > ParticleContainerType; - typedef vector< ParticleContainerType* > FiberContainerType; - vtkSmartPointer m_VtkCellArray; - vtkSmartPointer m_VtkPoints; + typedef itk::Image ItkFloatImageType; - typedef itk::Vector OdfVectorType; - typedef itk::Image ItkQBallImgType; - ItkQBallImgType::Pointer m_ItkQBallImage; - float m_FiberLength; - itk::Point m_LastPoint; + FiberBuilder(ParticleGrid* grid, ItkFloatImageType* image); + ~FiberBuilder(); - FiberBuilder(float *points, int numPoints, double spacing[], ItkQBallImgType::Pointer image); + vtkSmartPointer iterate(int minFiberLength); - ~FiberBuilder(); +protected: - vtkSmartPointer iterate(int minFiberLength); + void AddPoint(Particle *dp, vtkSmartPointer container); - void AddPoint(Particle *dp, vtkSmartPointer container); + void labelPredecessors(Particle *dp, vtkSmartPointer container); + void labelSuccessors(Particle *dp, vtkSmartPointer container); - void labelPredecessors(Particle *dp, vtkSmartPointer container); + itk::Point m_LastPoint; + float m_FiberLength; + ItkFloatImageType::Pointer m_Image; + ParticleGrid* m_Grid; + Particle* particles; + vtkSmartPointer m_VtkCellArray; + vtkSmartPointer m_VtkPoints; - void labelSuccessors(Particle *dp, vtkSmartPointer container); }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp index 266a54d546..194a357213 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp @@ -1,704 +1,646 @@ /*=================================================================== 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(MTRand* rgen, float *points,int numPoints, float *dimg, const int *dsz, double *voxsz, double cellsize) - : m_QBallImageData(dimg) - , datasz(dsz) - , enc(0) - , voxsize(voxsz) - , m_NumAttributes(0) +MetropolisHastingsSampler::MetropolisHastingsSampler(ParticleGrid* grid, MTRand* randGen) + : m_NumAttributes(0) , m_AcceptedProposals(0) { - mtrand = rgen; + mtrand = randGen; - m_ParticleGrid = new ParticleGrid(); - - width = dsz[1]*voxsize[0]; - height = dsz[2]*voxsize[1]; - depth = dsz[3]*voxsize[2]; - - fprintf(stderr,"Data dimensions (mm) : %f x %f x %f\n",width,height,depth); - fprintf(stderr,"Data dimensions (voxel) : %i x %i x %i\n",datasz[1],datasz[2],datasz[3]); - fprintf(stderr,"voxel size (mm) : %lf x %lf x %lf\n",voxsize[0],voxsize[1],voxsize[2]); - - float cellcnt_x = (int)((float)width/cellsize) +1; - float cellcnt_y = (int)((float)height/cellsize) +1; - float cellcnt_z = (int)((float)depth/cellsize) +1; - //int cell_capacity = 2048; - //int cell_capacity = 64; - int cell_capacity = 512; - - fprintf(stderr,"grid dimensions : %f x %f x %f\n",cellcnt_x,cellcnt_y,cellcnt_z); - fprintf(stderr,"grid cell size (mm) : %f^3\n",cellsize); - fprintf(stderr,"cell capacity : %i\n",cell_capacity); - fprintf(stderr,"#cells*cellcap : %.1f K\n",cell_capacity*cellcnt_x*cellcnt_y*cellcnt_z/1000); - - int minsize = 1000000; - int err = m_ParticleGrid->allocate(((numPoints>minsize)? (numPoints+100000) : minsize), cellcnt_x, cellcnt_y, cellcnt_z, cellsize, cell_capacity); - - if (err == -1) - { - fprintf(stderr,"MetropolisHastingsSamplerBase: out of Memory!\n"); - return; - } - - m_NumAttributes = 10; - for (int k = 0; k < numPoints; k++) - { - Particle *p = m_ParticleGrid->newParticle(vnl_vector_fixed(points[m_NumAttributes*k], points[m_NumAttributes*k+1],points[m_NumAttributes*k+2])); - if (p!=0) - { - p->N = vnl_vector_fixed(points[m_NumAttributes*k+3],points[m_NumAttributes*k+4],points[m_NumAttributes*k+5]); - p->cap = points[m_NumAttributes*k+6]; - p->len = points[m_NumAttributes*k+7]; - p->mID = (int) points[m_NumAttributes*k+8]; - p->pID = (int) points[m_NumAttributes*k+9]; - if (p->mID != -1) - m_ParticleGrid->m_NumConnections++; - if (p->pID != -1) - m_ParticleGrid->m_NumConnections++; - p->label = 0; - } - else - { - fprintf(stderr,"error: cannot allocate particle, con. indices will be wrong! \n"); - } - } - m_ParticleGrid->m_NumConnections /= 2; + m_ParticleGrid = grid; m_Iterations = 0; m_AcceptedProposals = 0; - externalEnergy = 0; internalEnergy = 0; } void MetropolisHastingsSampler::SetEnergyComputer(EnergyComputer *e) { enc = e; } void MetropolisHastingsSampler::Iterate(float* acceptance, unsigned long* numCon, unsigned long* numPart, bool *abort) { m_AcceptedProposals = 0; for (int it = 0; it < m_Iterations;it++) { if (*abort) break; IterateOneStep(); *numCon = m_ParticleGrid->m_NumConnections; *numPart = m_ParticleGrid->m_NumParticles; } *acceptance = (float)m_AcceptedProposals/m_Iterations; } void MetropolisHastingsSampler::WriteOutParticles(float *npoints) { for (int k = 0; k < m_ParticleGrid->m_NumParticles; k++) { Particle *p = &(m_ParticleGrid->m_Particles[k]); npoints[m_NumAttributes*k] = p->R[0]; npoints[m_NumAttributes*k+1] = p->R[1]; npoints[m_NumAttributes*k+2] = p->R[2]; npoints[m_NumAttributes*k+3] = p->N[0]; npoints[m_NumAttributes*k+4] = p->N[1]; npoints[m_NumAttributes*k+5] = p->N[2]; npoints[m_NumAttributes*k+6] = p->cap; npoints[m_NumAttributes*k+7] = p->len; - npoints[m_NumAttributes*k+8] = m_ParticleGrid->ID_2_index(p->mID); - npoints[m_NumAttributes*k+9] = m_ParticleGrid->ID_2_index(p->pID); + npoints[m_NumAttributes*k+8] = m_ParticleGrid->Id2Index(p->mID); + npoints[m_NumAttributes*k+9] = m_ParticleGrid->Id2Index(p->pID); } } void MetropolisHastingsSampler::SetParameters(float Temp, int numit, float plen, float curv_hardthres, float chempot_particle) { m_Iterations = numit; p_birth = 0.25; p_death = 0.05; p_shift = 0.15; p_shiftopt = 0.1; p_con = 0.45; p_cap = 0.0; m_ChempotParticle = chempot_particle; float sum = p_birth+p_death+p_shift+p_shiftopt+p_con; p_birth /= sum; p_death /= sum; p_shift /= sum; p_shiftopt /= sum; T_in = Temp; T_ex = 0.01; dens = exp(-chempot_particle/T_in); len_def = plen; len_sig = 0.0; cap_def = 1.0; cap_sig = 0.0; // shift proposal sigma_g = len_def/8.0; gamma_g = 1/(sigma_g*sigma_g*2); Z_g = pow(2*M_PI*sigma_g,3.0/2.0)*(M_PI*sigma_g/len_def); // conn proposal dthres = len_def; nthres = curv_hardthres; T_prop = 0.5; dthres *= dthres; stopprobability = exp(-1/T_prop); del_prob = 0.1; } void MetropolisHastingsSampler::SetTemperature(float temp) { T_in = temp; dens = exp(-m_ChempotParticle/T_in); } vnl_vector_fixed MetropolisHastingsSampler::distortn(float sigma, vnl_vector_fixed& vec) { vec[0] += sigma*mtrand->frandn(); vec[1] += sigma*mtrand->frandn(); vec[2] += sigma*mtrand->frandn(); } vnl_vector_fixed MetropolisHastingsSampler::rand_sphere() { vnl_vector_fixed vec; vec[0] += mtrand->frandn(); vec[1] += mtrand->frandn(); vec[2] += mtrand->frandn(); vec.normalize(); return vec; } //vnl_vector_fixed MetropolisHastingsSampler::rand(float w, float h, float d) //{ // vnl_vector_fixed vec; // vec[0] = mtrand->frand()*w; // vec[1] = mtrand->frand()*h; // vec[2] = mtrand->frand()*d; // return vec; //} void MetropolisHastingsSampler::IterateOneStep() { float randnum = mtrand->frand(); //randnum = 0; /////////////////////////////////////////////////////////////// //////// Birth Proposal /////////////////////////////////////////////////////////////// if (randnum < p_birth) { #ifdef TIMING tic(&birthproposal_time); birthstats.propose(); #endif vnl_vector_fixed R; enc->drawSpatPosition(R); //fprintf(stderr,"drawn: %f, %f, %f\n",R[0],R[1],R[2]); //R.setXYZ(20.5*3, 35.5*3, 1.5*3); vnl_vector_fixed N = rand_sphere(); //N.setXYZ(1,0,0); float cap = cap_def - cap_sig*mtrand->frand(); float len = len_def;// + len_sig*(mtrand->frand()-0.5); Particle prop; prop.R = R; prop.N = N; prop.cap = cap; prop.len = len; float prob = dens * p_death /((p_birth)*(m_ParticleGrid->m_NumParticles+1)); float ex_energy = enc->computeExternalEnergy(R,N,cap,len,0); float in_energy = enc->computeInternalEnergy(&prop); prob *= exp((in_energy/T_in+ex_energy/T_ex)) ; if (prob > 1 || mtrand->frand() < prob) { Particle *p = m_ParticleGrid->newParticle(R); if (p!=0) { p->R = R; p->N = N; p->cap = cap; p->len = len; #ifdef TIMING birthstats.accepted(); #endif m_AcceptedProposals++; } } #ifdef TIMING toc(&birthproposal_time); #endif } /////////////////////////////////////////////////////////////// //////// Death Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death) { if (m_ParticleGrid->m_NumParticles > 0) { #ifdef TIMING tic(&deathproposal_time); deathstats.propose(); #endif int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *dp = &(m_ParticleGrid->m_Particles[pnum]); if (dp->pID == -1 && dp->mID == -1) { float ex_energy = enc->computeExternalEnergy(dp->R,dp->N,dp->cap,dp->len,dp); float in_energy = enc->computeInternalEnergy(dp); float prob = m_ParticleGrid->m_NumParticles * (p_birth) /(dens*p_death); //*SpatProb(dp->R); prob *= exp(-(in_energy/T_in+ex_energy/T_ex)) ; if (prob > 1 || mtrand->frand() < prob) { m_ParticleGrid->RemoveParticle(pnum); #ifdef TIMING deathstats.accepted(); #endif m_AcceptedProposals++; } } #ifdef TIMING toc(&deathproposal_time); #endif } } /////////////////////////////////////////////////////////////// //////// Cap change Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death+p_cap) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); Particle prop_p = *p; prop_p.cap = cap_def - cap_sig*mtrand->frand(); float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,prop_p.cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); //float in_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) // - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float prob = exp(ex_energy/T_ex); // * SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand->frand() < prob) { p->cap = prop_p.cap; m_AcceptedProposals++; } } } /////////////////////////////////////////////////////////////// //////// Shift Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death+p_shift+p_cap) { float energy = 0; if (m_ParticleGrid->m_NumParticles > 0) { #ifdef TIMING tic(&shiftproposal_time); shiftstats.propose(); #endif int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); Particle prop_p = *p; distortn(sigma_g, prop_p.R); distortn(sigma_g/(2*p->len), prop_p.N); prop_p.N.normalize(); float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float in_energy = enc->computeInternalEnergy(&prop_p) - enc->computeInternalEnergy(p); float prob = exp(ex_energy/T_ex+in_energy/T_in); // * SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand->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; } #ifdef TIMING shiftstats.accepted(); #endif m_AcceptedProposals++; } #ifdef TIMING toc(&shiftproposal_time); #endif } } else if (randnum < p_birth+p_death+p_shift+p_shiftopt+p_cap) { float energy = 0; if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); bool no_proposal = false; Particle prop_p = *p; if (p->pID != -1 && p->mID != -1) { Particle *plus = m_ParticleGrid->m_AddressContainer[p->pID]; int ep_plus = (plus->pID == p->ID)? 1 : -1; Particle *minus = m_ParticleGrid->m_AddressContainer[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 *= 0.5; prop_p.N = plus->R - minus->R; prop_p.N.normalize(); } else if (p->pID != -1) { Particle *plus = m_ParticleGrid->m_AddressContainer[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.N = plus->N; } else if (p->mID != -1) { Particle *minus = m_ParticleGrid->m_AddressContainer[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.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))*gamma_g)/Z_g; float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float in_energy = enc->computeInternalEnergy(&prop_p) - enc->computeInternalEnergy(p); float prob = exp(ex_energy/T_ex+in_energy/T_in)*p_shift*p_rev/(p_shiftopt+p_shift*p_rev); //* SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand->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) { #ifdef TIMING tic(&connproposal_time); connstats.propose(); #endif int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); EndPoint P; P.p = p; P.ep = (mtrand->frand() > 0.5)? 1 : -1; RemoveAndSaveTrack(P); if (TrackBackup.proposal_probability != 0) { MakeTrackProposal(P); float prob = (TrackProposal.energy-TrackBackup.energy)/T_in ; // prob = exp(prob)*(TrackBackup.proposal_probability) // /(TrackProposal.proposal_probability); prob = exp(prob)*(TrackBackup.proposal_probability * pow(del_prob,TrackProposal.length)) /(TrackProposal.proposal_probability * pow(del_prob,TrackBackup.length)); if (mtrand->frand() < prob) { ImplementTrack(TrackProposal); #ifdef TIMING connstats.accepted(); #endif m_AcceptedProposals++; } else { ImplementTrack(TrackBackup); } } else ImplementTrack(TrackBackup); #ifdef TIMING toc(&connproposal_time); #endif } } } void MetropolisHastingsSampler::ImplementTrack(Track &T) { for (int k = 1; k < T.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; TrackBackup.track[cnt] = Current; EndPoint Next; for (;;) { Next.p = 0; if (Current.ep == 1) { if (Current.p->pID != -1) { Next.p = m_ParticleGrid->m_AddressContainer[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->m_AddressContainer[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 *= (simpsamp.probFor(Next)); if (Next.p == 0) // no successor -> break break; energy += enc->computeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); Current = Next; Current.ep *= -1; cnt++; TrackBackup.track[cnt] = Current; if (mtrand->rand() > del_prob) { break; } } TrackBackup.energy = energy; TrackBackup.proposal_probability = AccumProb; TrackBackup.length = cnt+1; } void MetropolisHastingsSampler::MakeTrackProposal(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; TrackProposal.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); // // no candidates anymore // if (simpsamp.isempty()) // break; int k = simpsamp.draw(mtrand->frand()); // stop tracking proposed if (k==0) break; EndPoint Next = simpsamp.objs[k]; float probability = simpsamp.probFor(k); // accumulate energy and proposal distribution energy += enc->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 TrackProposal.track[cnt++] = Current; } TrackProposal.energy = energy; TrackProposal.proposal_probability = AccumProb; TrackProposal.length = cnt; // clear labels for (int j = 0; j < TrackProposal.length;j++) TrackProposal.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) ); m_ParticleGrid->ComputeNeighbors(R); simpsamp.clear(); simpsamp.add(stopprobability,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(); if (dist < dthres) { dot = dot_product(p2->N,p->N) * ep; if (dot > nthres) { float en = enc->computeInternalEnergyConnection(p,ep,p2,-1); simpsamp.add(exp(en/T_prop),EndPoint(p2,-1)); } } } if (p2->pID == -1) { dist = (p2->R + p2->N * p2->len - R).squared_magnitude(); if (dist < dthres) { dot = dot_product(p2->N,p->N) * (-ep); if (dot > nthres) { float en = enc->computeInternalEnergyConnection(p,ep,p2,+1); simpsamp.add(exp(en/T_prop),EndPoint(p2,+1)); } } } } } } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h index 3df94e930a..af674b2e50 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h @@ -1,107 +1,113 @@ /*=================================================================== 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 _SAMPLER #define _SAMPLER +// MITK #include #include "mitkParticleGrid.h" #include "mitkEnergyComputer.h" #include "MersenneTwister.h" #include "mitkSimpSamp.h" + +// ITK +#include + +// MISC #include namespace mitk { class MitkDiffusionImaging_EXPORT MetropolisHastingsSampler { public: + typedef itk::Image< float, 3 > ItkFloatImageType; + ParticleGrid* m_ParticleGrid; - float* m_QBallImageData; const int* datasz; EnergyComputer* enc; int m_Iterations; float width; float height; float depth; - double *voxsize; int m_NumAttributes; int m_AcceptedProposals; float T_in ; float T_ex ; float dens; float p_birth; float p_death; float p_shift; float p_shiftopt; float p_cap; float p_con; float sigma_g; float gamma_g; float Z_g; float dthres; float nthres; float T_prop; float stopprobability; float del_prob; float len_def; float len_sig; float cap_def; float cap_sig; float externalEnergy; float internalEnergy; float m_ChempotParticle; MTRand* mtrand; Track TrackProposal, TrackBackup; SimpSamp simpsamp; - MetropolisHastingsSampler(MTRand* rgen, float *points,int numPoints, float *dimg, const int *dsz, double *voxsz, double cellsize); + MetropolisHastingsSampler(ParticleGrid* grid, MTRand* randGen); void WriteOutParticles(float *npoints); void SetEnergyComputer(EnergyComputer *e); void SetParameters(float Temp, int numit, float plen, float curv_hardthres, float chempot_particle); void SetTemperature(float temp); void Iterate(float* acceptance, unsigned long* numCon, unsigned long* numPart, bool *abort); void IterateOneStep(); void ImplementTrack(Track &T); void RemoveAndSaveTrack(EndPoint P); void MakeTrackProposal(EndPoint P); void ComputeEndPointProposalDistribution(EndPoint P); vnl_vector_fixed distortn(float sigma, vnl_vector_fixed& vec); vnl_vector_fixed rand_sphere(); // vnl_vector_fixed rand(float w, float h, float d); }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp index ce7afad1cb..1db8dcf8c0 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp @@ -1,430 +1,424 @@ /*=================================================================== 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 "mitkParticleGrid.h" #include #include using namespace mitk; -ParticleGrid::ParticleGrid() +ParticleGrid::ParticleGrid(ItkFloatImageType* image, float cellSize) { //// involving the container capacity = 0; m_Particles = 0; m_AddressContainer = 0; m_NumParticles = 0; m_NumConnections = 0; m_NumCellOverflows = 0; ////// involvin the grid nx = 0; ny = 0; nz = 0; csize = 0; gridsize = 0; grid = (Particle**) 0; occnt = (int*) 0; increase_step = 100000; m_Memory = 0; -} -float ParticleGrid::GetMemoryUsage() -{ - return m_Memory; -} + int width = image->GetLargestPossibleRegion().GetSize()[0]*image->GetSpacing()[0]; + int height = image->GetLargestPossibleRegion().GetSize()[1]*image->GetSpacing()[1]; + int depth = image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]; + + float cellcnt_x = (int)((float)width/cellSize) +1; + float cellcnt_y = (int)((float)height/cellSize) +1; + float cellcnt_z = (int)((float)depth/cellSize) +1; + int cell_capacity = 512; -int ParticleGrid::allocate(int _capacity, int _nx, int _ny, int _nz, float cellsize, int cellcapacity) -{ //// involving the container - capacity = _capacity; + capacity = 100000; m_Particles = (Particle*) malloc(sizeof(Particle)*capacity); m_AddressContainer = (Particle**) malloc(sizeof(Particle*)*capacity); if (m_Particles == 0 || m_AddressContainer == 0) { fprintf(stderr,"error: Out of Memory\n"); capacity = 0; - return -1; } else { fprintf(stderr,"Allocated Memory for %i particles \n",capacity); } - m_NumParticles = 0; int i; for (i = 0;i < capacity;i++) { - m_AddressContainer[i] = &(m_Particles[i]); // initialize pointer in LUT - m_Particles[i].ID = i; // initialize unique IDs + m_AddressContainer[i] = &(m_Particles[i]); // initialize pointer in LUT + m_Particles[i].ID = i; // initialize unique IDs } - ////// involvin the grid - nx = _nx; ny = _ny; nz = _nz; csize = cellcapacity; + // involving the grid + nx = cellcnt_x; ny = cellcnt_y; nz = cellcnt_z; csize = cell_capacity; gridsize = nx*ny*nz; m_Memory = (float)(sizeof(Particle*)*gridsize*csize)/1000000; grid = (Particle**) malloc(sizeof(Particle*)*gridsize*csize); occnt = (int*) malloc(sizeof(int)*gridsize); if (grid == 0 || occnt == 0) { fprintf(stderr,"error: Out of Memory\n"); capacity = 0; - return -1; } for (i = 0;i < gridsize;i++) occnt[i] = 0; - mulx = 1/cellsize; - muly = 1/cellsize; - mulz = 1/cellsize; - - return 1; + mulx = 1/cellSize; + muly = 1/cellSize; + mulz = 1/cellSize; } - - int ParticleGrid::reallocate() { int new_capacity = capacity + increase_step; Particle* new_particles = (Particle*) realloc(m_Particles,sizeof(Particle)*new_capacity); Particle** new_ID_2_address = (Particle**) realloc(m_AddressContainer,sizeof(Particle*)*new_capacity); if (new_particles == 0 || new_ID_2_address == 0) { fprintf(stderr,"ParticleGird:reallocate: out of memory!\n"); return -1; } fprintf(stderr," now %i particles are allocated \n",new_capacity); m_Memory = (float)(sizeof(Particle*)*new_capacity)/1000000; int i; for (i = 0; i < capacity; i++) { new_ID_2_address[i] = new_ID_2_address[i] - m_Particles + new_particles; // shift address } for (i = capacity; i < new_capacity; i++) { new_particles[i].ID = i; // initialize unique IDs new_ID_2_address[i] = &(new_particles[i]) ; // initliaze pointer in LUT } for (i = 0; i < nx*ny*nz*csize; i++) { grid[i] = grid[i] - m_Particles + new_particles; } m_Particles = new_particles; m_AddressContainer = new_ID_2_address; capacity = new_capacity; return 1; } ParticleGrid::~ParticleGrid() { if (m_Particles != 0) free(m_Particles); if (grid != 0) free(grid); if (occnt != 0) free(occnt); if (m_AddressContainer != 0) free(m_AddressContainer); } -int ParticleGrid::ID_2_index(int ID) +int ParticleGrid::Id2Index(int ID) { if (ID == -1) return -1; else return (m_AddressContainer[ID] - m_Particles); } Particle* ParticleGrid::newParticle(vnl_vector_fixed R) { /////// get free place in container; if (m_NumParticles >= capacity) { fprintf(stderr,"capacity overflow , reallocating ...\n"); if (reallocate() == -1) { fprintf(stderr,"out of Memory!!\n"); return 0; } } int xint = int(R[0]*mulx); if (xint < 0) { //fprintf(stderr,"error: out of grid\n"); return 0;} if (xint >= nx) { // fprintf(stderr,"error: out of grid\n"); return 0;} int yint = int(R[1]*muly); if (yint < 0) { //fprintf(stderr,"error: out of grid\n"); return 0;} if (yint >= ny) {// fprintf(stderr,"error: out of grid\n"); return 0;} int zint = int(R[2]*mulz); if (zint < 0) {// fprintf(stderr,"error: out of grid\n"); return 0;} if (zint >= nz) { //fprintf(stderr,"error: out of grid\n"); return 0;} int idx = xint + nx*(yint + ny*zint); if (occnt[idx] < csize) { Particle *p = &(m_Particles[m_NumParticles]); p->R = R; p->mID = -1; p->pID = -1; m_NumParticles++; p->gridindex = csize*idx + occnt[idx]; grid[p->gridindex] = p; occnt[idx]++; return p; } else { m_NumCellOverflows++; //fprintf(stderr,"error: cell overflow \n"); return 0; } } bool ParticleGrid::TryUpdateGrid(int k) { Particle* p = &(m_Particles[k]); /////// find new grid cell int xint = int(p->R[0]*mulx); if (xint < 0) { return false; } if (xint >= nx) { return false; } int yint = int(p->R[1]*muly); if (yint < 0) { return false; } if (yint >= ny) { return false; } int zint = int(p->R[2]*mulz); if (zint < 0) { return false; } if (zint >= nz) { return false; } int idx = xint + nx*(yint+ zint*ny); int cellidx = p->gridindex/csize; if (idx != cellidx) // cell has changed { if (occnt[idx] < csize) { // remove from old position in grid; int grdindex = p->gridindex; grid[grdindex] = grid[cellidx*csize + occnt[cellidx]-1]; grid[grdindex]->gridindex = grdindex; occnt[cellidx]--; // insert at new position in grid p->gridindex = idx*csize + occnt[idx]; grid[p->gridindex] = p; occnt[idx]++; return true; } else return false; } return true; } void ParticleGrid::RemoveParticle(int k) { Particle* p = &(m_Particles[k]); int grdindex = p->gridindex; int cellidx = grdindex/csize; int idx = grdindex%csize; // remove pending connections if (p->mID != -1) DestroyConnection(p,-1); if (p->pID != -1) DestroyConnection(p,+1); // remove from grid if (idx < occnt[cellidx]-1) { grid[grdindex] = grid[cellidx*csize + occnt[cellidx]-1]; grid[grdindex]->gridindex = grdindex; } occnt[cellidx]--; // remove from container if (kID; int move_ID = m_Particles[m_NumParticles-1].ID; *p = m_Particles[m_NumParticles-1]; // move very last particle to empty slot m_Particles[m_NumParticles-1].ID = todel_ID; // keep IDs unique grid[p->gridindex] = p; // keep gridindex consistent // permute address table m_AddressContainer[todel_ID] = &(m_Particles[m_NumParticles-1]); m_AddressContainer[move_ID] = p; } m_NumParticles--; } void ParticleGrid::ComputeNeighbors(vnl_vector_fixed &R) { float xfrac = R[0]*mulx; float yfrac = R[1]*muly; float zfrac = R[2]*mulz; int xint = int(xfrac); int yint = int(yfrac); int zint = int(zfrac); int dx = -1; if (xfrac-xint > 0.5) dx = 1; if (xint <= 0) { xint = 0; dx = 1; } if (xint >= nx-1) { xint = nx-1; dx = -1; } int dy = -1; if (yfrac-yint > 0.5) dy = 1; if (yint <= 0) {yint = 0; dy = 1; } if (yint >= ny-1) {yint = ny-1; dy = -1;} int dz = -1; if (zfrac-zint > 0.5) dz = 1; if (zint <= 0) {zint = 0; dz = 1; } if (zint >= nz-1) {zint = nz-1; dz = -1;} m_NeighbourTracker.cellidx[0] = xint + nx*(yint+zint*ny); m_NeighbourTracker.cellidx[1] = m_NeighbourTracker.cellidx[0] + dx; m_NeighbourTracker.cellidx[2] = m_NeighbourTracker.cellidx[1] + dy*nx; m_NeighbourTracker.cellidx[3] = m_NeighbourTracker.cellidx[2] - dx; m_NeighbourTracker.cellidx[4] = m_NeighbourTracker.cellidx[0] + dz*nx*ny; m_NeighbourTracker.cellidx[5] = m_NeighbourTracker.cellidx[4] + dx; m_NeighbourTracker.cellidx[6] = m_NeighbourTracker.cellidx[5] + dy*nx; m_NeighbourTracker.cellidx[7] = m_NeighbourTracker.cellidx[6] - dx; m_NeighbourTracker.cellidx_c[0] = csize*m_NeighbourTracker.cellidx[0]; m_NeighbourTracker.cellidx_c[1] = csize*m_NeighbourTracker.cellidx[1]; m_NeighbourTracker.cellidx_c[2] = csize*m_NeighbourTracker.cellidx[2]; m_NeighbourTracker.cellidx_c[3] = csize*m_NeighbourTracker.cellidx[3]; m_NeighbourTracker.cellidx_c[4] = csize*m_NeighbourTracker.cellidx[4]; m_NeighbourTracker.cellidx_c[5] = csize*m_NeighbourTracker.cellidx[5]; m_NeighbourTracker.cellidx_c[6] = csize*m_NeighbourTracker.cellidx[6]; m_NeighbourTracker.cellidx_c[7] = csize*m_NeighbourTracker.cellidx[7]; m_NeighbourTracker.cellcnt = 0; m_NeighbourTracker.pcnt = 0; } Particle* ParticleGrid::GetNextNeighbor() { if (m_NeighbourTracker.pcnt < occnt[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]]) { return grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt] + (m_NeighbourTracker.pcnt++)]; } else { for(;;) { m_NeighbourTracker.cellcnt++; if (m_NeighbourTracker.cellcnt >= 8) return 0; if (occnt[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]] > 0) break; } m_NeighbourTracker.pcnt = 1; return grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt]]; } } void ParticleGrid::CreateConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = P2->ID; else P1->pID = P2->ID; if (ep2 == -1) P2->mID = P1->ID; else P2->pID = P1->ID; m_NumConnections++; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = -1; else P1->pID = -1; if (ep2 == -1) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1) { Particle *P2 = 0; if (ep1 == 1) { P2 = m_AddressContainer[P1->pID]; P1->pID = -1; } else { P2 = m_AddressContainer[P1->mID]; P1->mID = -1; } if (P2->mID == P1->ID) { P2->mID = -1; } else { P2->pID = -1; } m_NumConnections--; } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h index 4c4b48b8bb..e69b74381d 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h @@ -1,118 +1,123 @@ /*=================================================================== 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 _PARTICLEGRID #define _PARTICLEGRID +// MITK #include "mitkParticle.h" #include "MitkDiffusionImagingExports.h" +// ITK +#include + namespace mitk { class MitkDiffusionImaging_EXPORT ParticleGrid { //////////////// Container public: + + typedef itk::Image< float, 3 > ItkFloatImageType; + Particle* m_Particles; // particles in linear array Particle** m_AddressContainer; int m_NumParticles; // number of particles int m_NumConnections; // number of connections int m_NumCellOverflows; // number of cell overflows - ParticleGrid(); + ParticleGrid(ItkFloatImageType* image, float cellSize); ~ParticleGrid(); - float GetMemoryUsage(); - int allocate(int _capacity, int _nx, int _ny, int _nz, float cellsize, int cellcapacity); int reallocate(); - int ID_2_index(int ID); + int Id2Index(int ID); Particle* newParticle(vnl_vector_fixed R); bool TryUpdateGrid(int k); void RemoveParticle(int k); void ComputeNeighbors(vnl_vector_fixed &R); Particle* GetNextNeighbor(); void CreateConnection(Particle *P1,int ep1, Particle *P2, int ep2); void DestroyConnection(Particle *P1,int ep1, Particle *P2, int ep2); void DestroyConnection(Particle *P1,int ep1); private: int capacity; // maximal number of particles int increase_step; Particle **grid; // the grid // grid size int nx; int ny; int nz; // scaling factor for grid float mulx; float muly; float mulz; int csize; // particle capacity of single cell in grid int *occnt; // occupation count of grid cells int gridsize; // total number of cells float m_Memory; struct NeighborTracker // to run over the neighbors { int cellidx[8]; int cellidx_c[8]; int cellcnt; int pcnt; } m_NeighbourTracker; }; class MitkDiffusionImaging_EXPORT Track { public: EndPoint track[1000]; float energy; float proposal_probability; int length; void clear() { length = 0; energy = 0; proposal_probability = 1; } bool isequal(Track &t) { for (int i = 0; i < length;i++) { if (track[i].p != t.track[i].p || track[i].ep != t.track[i].ep) return false; } return true; } }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index 16c92c3010..c4e94085b9 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,666 +1,380 @@ /*=================================================================== 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" #include #include +#include "GibbsTracking/MersenneTwister.h" #include "GibbsTracking/mitkFiberBuilder.h" +#include "GibbsTracking/mitkMetropolisHastingsSampler.h" +#include "GibbsTracking/mitkEnergyComputer.h" #pragma GCC visibility push(default) #include #pragma GCC visibility pop #include #include #include #include #include #include #include #include #include #include #include -struct LessDereference { - template - bool operator()(const T * lhs, const T * rhs) const { - return *lhs < *rhs; - } -}; - namespace itk{ -template< class TInputOdfImage, class TInputROIImage > -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -::GibbsTrackingFilter(): +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_ChempotParticle(0), 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_SubtractMean(true), m_BuildFibers(false), - m_Sampler(NULL), m_Steps(10), m_Memory(0), m_ProposalAcceptance(0), - m_GfaImage(NULL), m_CurvatureHardThreshold(0.7), m_Meanval_sq(0.0) { - //this->m_MeasurementFrame.set_identity(); - this->SetNumberOfRequiredInputs(2); //Filter needs a DWI image + a Mask Image -} -template< class TInputOdfImage, class TInputROIImage > -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -::~GibbsTrackingFilter(){ - if (m_Sampler!=NULL) - delete m_Sampler; } -//template< class TInputOdfImage, class TInputROIImage > -//void -//GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -//::ComputeFiberCorrelationOriginal(){ -// float bD = 15; - -// vnl_matrix_fixed bDir = -// *itk::PointShell >::DistributePointShell(); - -// const int N = QBALL_ODFSIZE; - -// vnl_matrix_fixed C = bDir.transpose()*bDir; -// vnl_matrix_fixed Q = C; -// 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; i 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; -// 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); - -// // OLD -// BESSEL_APPROXCOEFF[0] = 0,1982; -// BESSEL_APPROXCOEFF[1] = 0.3415; -// BESSEL_APPROXCOEFF[2] = -0.9515; -// BESSEL_APPROXCOEFF[3] = 1.3423; -//} - -template< class TInputOdfImage, class TInputROIImage > -void -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -::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; -} - -// build fibers from tracking result -template< class TInputOdfImage, class TInputROIImage > -void -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -::BuildFibers(float* points, int numPoints) +template< class ItkQBallImageType > +GibbsTrackingFilter< ItkQBallImageType >::~GibbsTrackingFilter() { - double spacing[3]; - spacing[0] = m_ItkQBallImage->GetSpacing().GetElement(0); - spacing[1] = m_ItkQBallImage->GetSpacing().GetElement(1); - spacing[2] = m_ItkQBallImage->GetSpacing().GetElement(2); - - m_FiberPolyData = FiberPolyDataType::New(); - // initialize array of particles - FiberBuilder fiberBuilder(points, numPoints, spacing, m_ItkQBallImage); - // label the particles according to fiber affiliation and return polydata - m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); - m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); - - MITK_INFO << "itkGibbsTrackingFilter: " << m_NumAcceptedFibers << " accepted"; } // fill output fiber bundle datastructure -template< class TInputOdfImage, class TInputROIImage > -typename GibbsTrackingFilter< TInputOdfImage, TInputROIImage >::FiberPolyDataType -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -::GetFiberBundle() +template< class ItkQBallImageType > +typename GibbsTrackingFilter< ItkQBallImageType >::FiberPolyDataType GibbsTrackingFilter< ItkQBallImageType >::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } return m_FiberPolyData; } -// get memory allocated for particle grid -template< class TInputOdfImage, class TInputROIImage > -float -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -::GetMemoryUsage() -{ - if (m_Sampler!=NULL) - return m_Sampler->m_ParticleGrid->GetMemoryUsage(); - return 0; -} - -template< class TInputOdfImage, class TInputROIImage > +template< class ItkQBallImageType > bool -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > +GibbsTrackingFilter< ItkQBallImageType > ::EstimateParticleWeight() { MITK_INFO << "itkGibbsTrackingFilter: estimating particle weight"; typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); - gfaFilter->SetInput(m_ItkQBallImage); + gfaFilter->SetInput(m_QBallImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); - m_GfaImage = gfaFilter->GetOutput(); + ItkFloatImageType::Pointer gfaImage = gfaFilter->GetOutput(); float samplingStart = 1.0; float samplingStop = 0.66; - // copy GFA image (original should not be changed) - typedef itk::ImageDuplicator< GfaImageType > DuplicateFilterType; - DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); - duplicator->SetInputImage( m_GfaImage ); - duplicator->Update(); - m_GfaImage = duplicator->GetOutput(); + // GFA iterator + typedef ImageRegionIterator< ItkFloatImageType > GfaIteratorType; + GfaIteratorType gfaIt(gfaImage, gfaImage->GetLargestPossibleRegion() ); - //// GFA iterator //// - typedef ImageRegionIterator< GfaImageType > GfaIteratorType; - GfaIteratorType gfaIt(m_GfaImage, m_GfaImage->GetLargestPossibleRegion() ); + // Mask iterator + typedef ImageRegionConstIterator< ItkFloatImageType > MaskIteratorType; + MaskIteratorType mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); - //// Mask iterator //// - typedef ImageRegionConstIterator< MaskImageType > MaskIteratorType; - MaskIteratorType maskIt(m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); - - // set unmasked region of gfa image to 0 - gfaIt.GoToBegin(); - maskIt.GoToBegin(); - while( !gfaIt.IsAtEnd() ) - { - if(maskIt.Get()<=0) - gfaIt.Set(0); - ++gfaIt; - ++maskIt; - } - - // rescale gfa image to [0,1] - typedef itk::RescaleIntensityImageFilter< GfaImageType, GfaImageType > RescaleFilterType; - RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New(); - rescaleFilter->SetInput( m_GfaImage ); - rescaleFilter->SetOutputMaximum( samplingStart ); - rescaleFilter->SetOutputMinimum( 0 ); - rescaleFilter->Update(); - m_GfaImage = rescaleFilter->GetOutput(); - gfaIt = GfaIteratorType(m_GfaImage, m_GfaImage->GetLargestPossibleRegion() ); - - //// Input iterator //// - typedef ImageRegionConstIterator< InputQBallImageType > InputIteratorType; - InputIteratorType git(m_ItkQBallImage, m_ItkQBallImage->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) { - git.GoToBegin(); + it.GoToBegin(); + mit.GoToBegin(); gfaIt.GoToBegin(); + while( !gfaIt.IsAtEnd() ) { - if(gfaIt.Get()>thr) + if(gfaIt.Get()>thr && mit.Get()>0) { - itk::OrientationDistributionFunction odf(git.Get().GetDataPointer()); + itk::OrientationDistributionFunction odf(it.Get().GetDataPointer()); upper += odf.GetMaxValue()-odf.GetMeanValue(); - ++count; } + ++it; + ++mit; ++gfaIt; - ++git; } } + if (count>0) upper /= count; else return false; m_ParticleWeight = upper/6; return true; } // perform global tracking -template< class TInputOdfImage, class TInputROIImage > -void -GibbsTrackingFilter< TInputOdfImage, TInputROIImage > -::GenerateData(){ +template< class ItkQBallImageType > +void GibbsTrackingFilter< ItkQBallImageType >::GenerateData() +{ - // input qball image - m_ItkQBallImage = dynamic_cast(this->GetInput(0)); m_NumAcceptedFibers = 0; - // approximationscoeffizienten der - // teilchenkorrelationen im orientierungsraum - // 4er vektor - //ComputeFiberCorrelationOriginal(); - ComputeFiberCorrelation(); - // image sizes and spacing - int qBallImageSize[4] = {QBALL_ODFSIZE, - m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(0), - m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(1), - m_ItkQBallImage->GetLargestPossibleRegion().GetSize().GetElement(2)}; - double qBallImageSpacing[3] = {m_ItkQBallImage->GetSpacing().GetElement(0),m_ItkQBallImage->GetSpacing().GetElement(1),m_ItkQBallImage->GetSpacing().GetElement(2)}; + 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 (qBallImageSize[1]<3 || qBallImageSize[2]<3 || qBallImageSize[3]<3) + 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_ItkQBallImage->GetDirection().GetVnlMatrix(); + 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 image buffer - int bufferSize = qBallImageSize[0]*qBallImageSize[1]*qBallImageSize[2]*qBallImageSize[3]; - float* qBallImageBuffer = (float*) m_ItkQBallImage->GetBufferPointer(); - float* workingQballImage = new float[bufferSize]; - for (int i=0; i DuplicateFilterType; + typename DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); + duplicator->SetInputImage( m_QBallImage ); + duplicator->Update(); + typename ItkQBallImageType::Pointer qballImage = duplicator->GetOutput(); // perform mean subtraction on odfs - if (m_SubtractMean) + typedef ImageRegionIterator< ItkQBallImageType > InputIteratorType; + InputIteratorType it(qballImage, qballImage->GetLargestPossibleRegion() ); + it.GoToBegin(); + while (!it.IsAtEnd()) { - float sum = 0; - for (int i=0; i0 && i%qBallImageSize[0] == 0 && i>0) - { - sum /= qBallImageSize[0]; - for (int j=i-qBallImageSize[0]; j 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 { mask = 0; - maskImageSize[0] = qBallImageSize[1]; - maskImageSize[1] = qBallImageSize[2]; - maskImageSize[2] = qBallImageSize[3]; + maskImageSize[0] = imgSize[1]; + maskImageSize[1] = imgSize[2]; + maskImageSize[2] = imgSize[3]; } - int mask_oversamp_mult = maskImageSize[0]/qBallImageSize[1]; + int mask_oversamp_mult = maskImageSize[0]/imgSize[1]; // load lookuptable 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(); + 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; + 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(); + 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; + MITK_INFO << "itkGibbsTrackingFilter: unable to open indices file"; + m_AbortTracking = true; } // initialize sphere interpolator with lookuptables SphereInterpolator *sinterp = new SphereInterpolator(coords, ind, QBALL_ODFSIZE, 301, 0.5); // get paramters float minSpacing; - if(qBallImageSpacing[0]m_NumIt) { MITK_INFO << "itkGibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } if (m_CurvatureHardThreshold < mitk::eps) - m_CurvatureHardThreshold = 0; + m_CurvatureHardThreshold = 0; unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); MTRand randGen(1); srand(1); - // setup metropolis hastings sampler MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; - if (m_Sampler!=NULL) - delete m_Sampler; - m_Sampler = new MetropolisHastingsSampler(&randGen, NULL, 0, workingQballImage, qBallImageSize, qBallImageSpacing, cellsize); - // setup energy computer - MITK_INFO << "itkGibbsTrackingFilter: setting up Energy-computer"; - EnergyComputer encomp(&randGen, workingQballImage,qBallImageSize,qBallImageSpacing,sinterp,m_Sampler->m_ParticleGrid,mask,mask_oversamp_mult, directionMatrix); - encomp.setParameters(m_ParticleWeight,m_ParticleWidth,m_ChempotConnection*m_ParticleLength*m_ParticleLength,m_ParticleLength,m_CurvatureHardThreshold,m_InexBalance,m_Chempot2, m_Meanval_sq); - m_Sampler->SetEnergyComputer(&encomp); - m_Sampler->SetParameters(m_TempStart,singleIts,m_ParticleLength,m_CurvatureHardThreshold,m_ChempotParticle); + ParticleGrid* particleGrid = new ParticleGrid(m_MaskImage, 2*m_ParticleLength); + MetropolisHastingsSampler* sampler = new MetropolisHastingsSampler(particleGrid, &randGen); + EnergyComputer encomp(&randGen, qballImage, imgSize,spacing,sinterp,particleGrid,mask,mask_oversamp_mult, directionMatrix); + encomp.setParameters(m_ParticleWeight,m_ParticleWidth,m_ChempotConnection*m_ParticleLength*m_ParticleLength,m_ParticleLength,m_CurvatureHardThreshold,m_InexBalance,m_Chempot2, m_Meanval_sq); - MITK_INFO << "itkGibbsTrackingFilter: Iterations: " << m_NumIt; - MITK_INFO << "itkGibbsTrackingFilter: steps: " << m_Steps; - MITK_INFO << "itkGibbsTrackingFilter: Particle weight: " << m_ParticleWeight; - MITK_INFO << "itkGibbsTrackingFilter: Particle width: " << m_ParticleWidth; - MITK_INFO << "itkGibbsTrackingFilter: Particle length: " << m_ParticleLength; - MITK_INFO << "itkGibbsTrackingFilter: Min. fiber length: " << m_ParticleLength; - MITK_INFO << "itkGibbsTrackingFilter: Start temperature: " << m_TempStart; - MITK_INFO << "itkGibbsTrackingFilter: End temperature: " << m_TempEnd; - MITK_INFO << "itkGibbsTrackingFilter: Energy balance: " << m_InexBalance; - MITK_INFO << "itkGibbsTrackingFilter: Curvature threshold: " << m_CurvatureHardThreshold; + sampler->SetEnergyComputer(&encomp); + sampler->SetParameters(m_TempStart,singleIts,m_ParticleLength,m_CurvatureHardThreshold,m_ChempotParticle); // main loop for( int step = 0; step < m_Steps; step++ ) { - if (m_AbortTracking) - break; - - m_CurrentStep = step+1; - float temperature = m_TempStart * exp(alpha*(((1.0)*step)/((1.0)*m_Steps))); - - m_Sampler->SetTemperature(temperature); - m_Sampler->Iterate(&m_ProposalAcceptance, &m_NumConnections, &m_NumParticles, &m_AbortTracking); - - 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)step/m_Steps << "%"; - - if (m_BuildFibers) - { - int numPoints = m_Sampler->m_ParticleGrid->m_NumParticles; - float* points = new float[numPoints*m_Sampler->m_NumAttributes]; - m_Sampler->WriteOutParticles(points); - BuildFibers(points, numPoints); - delete points; - m_BuildFibers = false; - } + if (m_AbortTracking) + break; + + m_CurrentStep = step+1; + float temperature = m_TempStart * exp(alpha*(((1.0)*step)/((1.0)*m_Steps))); + + sampler->SetTemperature(temperature); + sampler->Iterate(&m_ProposalAcceptance, &m_NumConnections, &m_NumParticles, &m_AbortTracking); + + 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)step/m_Steps << "%"; + + if (m_BuildFibers) + { + FiberBuilder fiberBuilder(particleGrid, m_MaskImage); + m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); + m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); + m_BuildFibers = false; + } } - int numPoints = m_Sampler->m_ParticleGrid->m_NumParticles; - float* points = new float[numPoints*m_Sampler->m_NumAttributes]; - m_Sampler->WriteOutParticles(points); - BuildFibers(points, numPoints); - delete points; + + FiberBuilder fiberBuilder(particleGrid, m_MaskImage); + m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); + m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); delete sinterp; delete coords; delete ind; - delete workingQballImage; + delete sampler; + delete particleGrid; m_AbortTracking = true; m_BuildFibers = false; MITK_INFO << "itkGibbsTrackingFilter: done generate data"; } } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h index 5919f6e018..857abc8e50 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h @@ -1,198 +1,153 @@ /*=================================================================== 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 #include "itkProcessObject.h" #include "itkVectorContainer.h" #include "itkImage.h" -#include "GibbsTracking/MersenneTwister.h" -#include "GibbsTracking/mitkMetropolisHastingsSampler.h" -#include "GibbsTracking/mitkEnergyComputer.h" - #include #include #include #include #include #include #include namespace itk{ - template< class TInputQBallImage, class TInputROIImage > - class GibbsTrackingFilter : - public ProcessObject{ - public: +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 ) - /** Types for the DWI Input Image **/ - typedef TInputQBallImage InputQBallImageType; - - /** Types for the Mask Image **/ - typedef TInputROIImage MaskImageType; - typedef typename MaskImageType::Pointer MaskImageTypePointer; - - typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; - - typedef Image< float, 3 > GfaImageType; - typedef typename GfaImageType::Pointer GfaImageTypePointer; + 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 ) 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( ChempotParticle, float ) itkGetMacro( ChempotParticle, float ) itkSetMacro( InexBalance, float ) itkGetMacro( InexBalance, float ) itkSetMacro( Chempot2, float ) itkGetMacro( Chempot2, float ) itkSetMacro( FiberLength, int ) itkGetMacro( FiberLength, int ) itkSetMacro( AbortTracking, bool ) itkGetMacro( AbortTracking, bool ) itkSetMacro( CurrentStep, unsigned long ) itkGetMacro( CurrentStep, unsigned long ) - itkSetMacro( SubtractMean, bool) - itkGetMacro( SubtractMean, bool) - itkSetMacro( CurvatureHardThreshold, float) itkGetMacro( CurvatureHardThreshold, float) - /** Set/Get the Odf Input Image **/ - itkSetInputMacro(OdfImage, InputQBallImageType, 0) - itkGetInputMacro(OdfImage, InputQBallImageType, 0) - - /** Set/Get the Input mask image **/ - itkSetMacro(MaskImage, MaskImageTypePointer) - itkGetMacro(MaskImage, MaskImageTypePointer) - - itkSetMacro(GfaImage, GfaImageTypePointer) - itkGetMacro(GfaImage, GfaImageTypePointer) - itkGetMacro(NumParticles, unsigned long) itkGetMacro(NumConnections, unsigned long) itkGetMacro(NumAcceptedFibers, int) itkGetMacro(ProposalAcceptance, float) itkGetMacro(Steps, unsigned int) - /** Entry Point For the Algorithm: Is invoked when Update() is called - either directly or through itk pipeline propagation - **/ - void GenerateData(); + // input data + itkSetMacro(QBallImage, typename ItkQBallImageType::Pointer) + itkSetMacro(MaskImage, ItkFloatImageType::Pointer) - /** override the Process Object Update because we don't have a - dataobject as an outpgnome themeut. We can change this later by wrapping the - tractcontainer in a dataobject decorator and letting the Superclass - know about it. - **/ - struct StochasticTractGenerationCallbackStruct{ - Pointer Filter; - }; + void GenerateData(); virtual void Update(){ - this->GenerateData(); + this->GenerateData(); } FiberPolyDataType GetFiberBundle(); - float GetMemoryUsage(); - bool EstimateParticleWeight(); - protected: +protected: GibbsTrackingFilter(); virtual ~GibbsTrackingFilter(); - - void ComputeFiberCorrelation(); - void ComputeFiberCorrelationOriginal(); - - void BuildFibers(float* points, int numPoints); + bool EstimateParticleWeight(); // Input Images - typename InputQBallImageType::Pointer m_ItkQBallImage; - typename MaskImageType::Pointer m_MaskImage; - typename GfaImageType::Pointer m_GfaImage; + typename ItkQBallImageType::Pointer m_QBallImage; + typename ItkFloatImageType::Pointer m_MaskImage; // 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_ChempotParticle;// unbenutzt (immer null, wenn groesser dann insgesamt weniger teilchen) - float m_InexBalance; // gewichtung zwischen den lambdas - // -5 ... 5 -> nur intern ... nur extern,default 0 - float m_Chempot2; // typischerweise 0, - // korrektur fuer das geschaetzte integral + 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_ChempotParticle; // unbenutzt (immer null, wenn groesser dann insgesamt weniger teilchen) + 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; - bool m_SubtractMean; int m_NumAcceptedFibers; volatile bool m_BuildFibers; unsigned int m_Steps; float m_Memory; float m_ProposalAcceptance; float m_CurvatureHardThreshold; float m_Meanval_sq; - MetropolisHastingsSampler* m_Sampler; FiberPolyDataType m_FiberPolyData; unsigned long m_NumParticles; unsigned long m_NumConnections; - }; +}; } #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 b1b5eb8239..a0590f3e1b 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingView.cpp @@ -1,770 +1,769 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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->SetInput0(m_View->m_ItkQBallImage.GetPointer()); +// 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->SetSubtractMean(m_View->m_Controls->m_MeanSubtractionCheckbox->isChecked()); 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 (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+"°"); } 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 74776790f4..d7019ced8b 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,165 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) -Copyright (c) German Cancer Research Center, +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 +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 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 +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 GibbsTrackingFilterType; + 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; mitk::QBallImage::Pointer m_QBallImage; ItkQBallImgType::Pointer m_ItkQBallImage; // data nodes mitk::DataNode::Pointer m_QBallImageNode; 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 de4d981ce9..2558f22cca 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkGibbsTrackingViewControls.ui @@ -1,1045 +1,1013 @@ QmitkGibbsTrackingViewControls 0 0 463 1011 0 0 0 0 QmitkTemplate 0 9 3 9 3 Data Q-Ball 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 - - + + - auto + 0 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter - - + + - auto = 0.5 * min. spacing; sigma + + + + + + + + + + Particle Width: + + + + + + + Only fibers longer than specified are accepted. - 100 + 500 1 + + 40 + Qt::Horizontal QSlider::NoTicks - - + + - auto + Curvature Threshold: + + + + + + + 0.001 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter - - + + - automatic estimation from gfa map and q-ball data. - - - 0 - - - 1000 - - - 1 + - - 0 + + - - Qt::Horizontal + + - - true + + 40mm - - QSlider::NoTicks + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 0.1 Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter - - - - 1 - - - 100 - - - 1 - - - 10 - - - Qt::Horizontal - - - false + + + + - - false + + - - QSlider::NoTicks + + - - - - - 0.001 + auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter - - - - 1 - - - 99 - - - 1 - - - 10 - - - Qt::Horizontal - - - QSlider::NoTicks - - - - - + + - 0 - - - Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + Balance In/Ex Energy: IE Bias < 0 < EE Bias -50 50 1 Qt::Horizontal QSlider::NoTicks - - + + - - - - - - - - - - 40mm - - - Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter + automatic estimation from gfa map and q-ball data. - - - - - - Only fibers longer than specified are accepted. + + 0 - 500 + 1000 1 - 40 + 0 Qt::Horizontal + + true + QSlider::NoTicks - - + + - Particle Length: + Min. Fiber Length: - - + + - Particle Width: + Particle Length: + + + + + + + 1 + + + 99 + + + 1 + + + 10 + + + Qt::Horizontal + + + QSlider::NoTicks Particle Weight: Start Temperature: - - - - End Temperature: - - - - - + + - + Allow only fiber curvature values smaller than the selected threshold. - - + + 180 - - + + 1 - - Balance In/Ex Energy: + + 45 + + + Qt::Horizontal + + + QSlider::NoTicks - - + + - + auto = 1.5 * min. spacing; l - - + + 100 - - + + 1 - - Min. Fiber Length: + + Qt::Horizontal + + + QSlider::NoTicks - - - - true - + + - Use mean subtracted ODFs (recommended). + auto = 0.5 * min. spacing; sigma - - Subtract ODF Mean + + 100 - - true + + 1 - - - - Qt::Horizontal - - QSizePolicy::Fixed - - - - 60 - 20 - + + QSlider::NoTicks - + - - + + - Curvature Threshold: + 45° + + + Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter - - + + - 45° + auto Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter - - - - Allow only fiber curvature values smaller than the selected threshold. - - - 180 - - - 1 - - - 45 - - - Qt::Horizontal - - - QSlider::NoTicks + + + + End Temperature: - - - - auto = 1.5 * min. spacing; l + + + + 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