diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp index d2eece6761..7c8cbed799 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp @@ -1,452 +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, 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); 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); 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); 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); 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); 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); 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); 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); 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 EnergyComputer::computeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len, Particle *dp) { float m = SpatProb(R); if (m == 0) return -INFINITY; float Dn = evaluateODF(R,N,len); float Sn = 0; float Pn = 0; m_ParticleGrid->ComputeNeighbors(R); for (;;) { Particle *p = m_ParticleGrid->GetNextNeighbor(); if (p == 0) break; if (dp != p) { float dot = fabs(dot_product(N,p->N)); float bw = mbesseli0(dot); float dpos = (p->R-R).squared_magnitude(); float w = mexp(dpos*gamma_s); - Sn += w*(bw+chempot2)*p->cap ; + Sn += w*(bw+chempot2); w = mexp(dpos*gamma_reg_s); Pn += w*bw; } } float energy = 0; - energy += (2*(Dn/particle_weight-Sn) - (mbesseli0(1.0)+chempot2)*cap)*cap; + energy += 2*(Dn/particle_weight-Sn) - (mbesseli0(1.0)+chempot2); return energy*ex_strength; } float EnergyComputer::computeInternalEnergy(Particle *dp) { float energy = eigen_energy; if (dp->pID != -1) energy += computeInternalEnergyConnection(dp,+1); if (dp->mID != -1) energy += computeInternalEnergyConnection(dp,-1); return energy; } float EnergyComputer::computeInternalEnergyConnection(Particle *p1,int ep1) { Particle *p2 = 0; int ep2; if (ep1 == 1) p2 = m_ParticleGrid->GetParticle(p1->pID); else p2 = m_ParticleGrid->GetParticle(p1->mID); if (p2->mID == p1->ID) ep2 = -1; else if (p2->pID == p1->ID) ep2 = 1; else fprintf(stderr,"EnergyComputer_connec: Connections are inconsistent!\n"); if (p2 == 0) fprintf(stderr,"bug2"); return computeInternalEnergyConnection(p1,ep1,p2,ep2); } float EnergyComputer::computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) { if ((dot_product(p1->N,p2->N))*ep1*ep2 > -curv_hard) return -INFINITY; vnl_vector_fixed R1 = p1->R + (p1->N * (p1->len * ep1)); vnl_vector_fixed R2 = p2->R + (p2->N * (p2->len * ep2)); if ((R1-R2).squared_magnitude() > particle_length_sq) return -INFINITY; vnl_vector_fixed R = (p2->R + p1->R); R *= 0.5; if (SpatProb(R) == 0) return -INFINITY; float norm1 = (R1-R).squared_magnitude(); float norm2 = (R2-R).squared_magnitude(); float energy = (eigencon_energy-norm1-norm2)*in_strength; return energy; } float EnergyComputer::mbesseli0(float x) { // BESSEL_APPROXCOEFF[0] = -0.1714; // BESSEL_APPROXCOEFF[1] = 0.5332; // BESSEL_APPROXCOEFF[2] = -1.4889; // BESSEL_APPROXCOEFF[3] = 2.0389; float y = x*x; float erg = -0.1714; erg += y*0.5332; erg += y*y*-1.4889; erg += y*y*y*2.0389; return erg; } float EnergyComputer::mexp(float x) { return((x>=7.0) ? 0 : ((x>=5.0) ? (-0.0029*x+0.0213) : ((x>=3.0) ? (-0.0215*x+0.1144) : ((x>=2.0) ? (-0.0855*x+0.3064) : ((x>=1.0) ? (-0.2325*x+0.6004) : ((x>=0.5) ? (-0.4773*x+0.8452) : ((x>=0.0) ? (-0.7869*x+1.0000) : 1 ))))))); // return exp(-x); } //ComputeFiberCorrelation() //{ // float bD = 15; // vnl_matrix_fixed bDir = // *itk::PointShell >::DistributePointShell(); // const int N = QBALL_ODFSIZE; // vnl_matrix_fixed temp = bDir.transpose(); // vnl_matrix_fixed C = temp*bDir; // vnl_matrix_fixed Q = C; // vnl_vector_fixed mean; // for(int i=0; i repMean; // for (int i=0; i P = Q*Q; // std::vector pointer; // pointer.reserve(N*N); // double * start = C.data_block(); // double * end = start + N*N; // for (double * iter = start; iter != end; ++iter) // { // pointer.push_back(iter); // } // std::sort(pointer.begin(), pointer.end(), LessDereference()); // vnl_vector_fixed alpha; // vnl_vector_fixed beta; // for (int i=0; im_Meanval_sq = (sum*sum)/N; // vnl_vector_fixed alpha_0; // vnl_vector_fixed alpha_2; // vnl_vector_fixed alpha_4; // vnl_vector_fixed alpha_6; // for(int i=0; i T; // T.set_column(0,alpha_0); // T.set_column(1,alpha_2); // T.set_column(2,alpha_4); // T.set_column(3,alpha_6); // vnl_vector_fixed coeff = vnl_matrix_inverse(T).pinverse()*beta; // MITK_INFO << "itkGibbsTrackingFilter: Bessel oefficients: " << coeff; // BESSEL_APPROXCOEFF = new float[4]; // BESSEL_APPROXCOEFF[0] = coeff(0); // BESSEL_APPROXCOEFF[1] = coeff(1); // BESSEL_APPROXCOEFF[2] = coeff(2); // BESSEL_APPROXCOEFF[3] = coeff(3); // BESSEL_APPROXCOEFF[0] = -0.1714; // BESSEL_APPROXCOEFF[1] = 0.5332; // BESSEL_APPROXCOEFF[2] = -1.4889; // BESSEL_APPROXCOEFF[3] = 2.0389; //} diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h index 4ccb73684a..8ef31620c9 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h @@ -1,102 +1,100 @@ /*=================================================================== 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 +#include +#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; 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, 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 computeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, float &len, Particle *dp); float computeInternalEnergy(Particle *dp); float computeInternalEnergyConnection(Particle *p1,int ep1); float computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2); float mbesseli0(float x); float mexp(float x); }; #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h index 5a16e52623..69af3bb6ba 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h @@ -1,69 +1,67 @@ /*=================================================================== 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 - +// MITK #include -#include - -#include "mitkParticleGrid.h" +#include +// VTK #include #include #include #include #include #include -#include #include using namespace std; namespace mitk { class MitkDiffusionImaging_EXPORT FiberBuilder { public: typedef itk::Image ItkFloatImageType; FiberBuilder(ParticleGrid* grid, ItkFloatImageType* image); ~FiberBuilder(); vtkSmartPointer iterate(int minFiberLength); protected: void AddPoint(Particle *dp, vtkSmartPointer container); void labelPredecessors(Particle *dp, vtkSmartPointer container); void labelSuccessors(Particle *dp, vtkSmartPointer container); itk::Point m_LastPoint; float m_FiberLength; ItkFloatImageType::Pointer m_Image; ParticleGrid* m_Grid; vtkSmartPointer m_VtkCellArray; vtkSmartPointer m_VtkPoints; }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp index 5823e730cf..dbd159e376 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp @@ -1,628 +1,580 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkMetropolisHastingsSampler.h" using namespace mitk; MetropolisHastingsSampler::MetropolisHastingsSampler(ParticleGrid* grid, MTRand* randGen) - : m_NumAttributes(0) - , m_AcceptedProposals(0) + : m_AcceptedProposals(0) { mtrand = randGen; - 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::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_BirthProb = 0.25; + m_DeathProb = 0.05; + m_ShiftProb = 0.15; + m_OptShiftProb = 0.1; + m_ConnectionProb = 0.45; 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; + float sum = m_BirthProb+m_DeathProb+m_ShiftProb+m_OptShiftProb+m_ConnectionProb; + m_BirthProb /= sum; m_DeathProb /= sum; m_ShiftProb /= sum; m_OptShiftProb /= sum; - T_in = Temp; - T_ex = 0.01; - dens = exp(-chempot_particle/T_in); + m_InTemp = Temp; + m_ExTemp = 0.01; + m_Density = exp(-chempot_particle/m_InTemp); 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); + m_InTemp = temp; + m_Density = exp(-m_ChempotParticle/m_InTemp); } 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) + if (randnum < m_BirthProb) { #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 prob = m_Density * m_DeathProb /((m_BirthProb)*(m_ParticleGrid->m_NumParticles+1)); - float ex_energy = enc->computeExternalEnergy(R,N,cap,len,0); + float ex_energy = enc->computeExternalEnergy(R,N,len,0); float in_energy = enc->computeInternalEnergy(&prop); - prob *= exp((in_energy/T_in+ex_energy/T_ex)) ; + prob *= exp((in_energy/m_InTemp+ex_energy/m_ExTemp)) ; 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) + else if (randnum < m_BirthProb+m_DeathProb) { 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->GetParticle(pnum); if (dp->pID == -1 && dp->mID == -1) { - float ex_energy = enc->computeExternalEnergy(dp->R,dp->N,dp->cap,dp->len,dp); + float ex_energy = enc->computeExternalEnergy(dp->R,dp->N,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)) ; + float prob = m_ParticleGrid->m_NumParticles * (m_BirthProb) /(m_Density*m_DeathProb); //*SpatProb(dp->R); + prob *= exp(-(in_energy/m_InTemp+ex_energy/m_ExTemp)) ; if (prob > 1 || 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->GetParticle(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) + else if (randnum < m_BirthProb+m_DeathProb+m_ShiftProb) { 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->GetParticle(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 ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->len,p) + - enc->computeExternalEnergy(p->R,p->N,p->len,p); float in_energy = enc->computeInternalEnergy(&prop_p) - enc->computeInternalEnergy(p); - float prob = exp(ex_energy/T_ex+in_energy/T_in); + float prob = exp(ex_energy/m_ExTemp+in_energy/m_InTemp); // * 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) + else if (randnum < m_BirthProb+m_DeathProb+m_ShiftProb+m_OptShiftProb) { float energy = 0; if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = m_ParticleGrid->GetParticle(pnum); bool no_proposal = false; Particle prop_p = *p; if (p->pID != -1 && p->mID != -1) { Particle *plus = m_ParticleGrid->GetParticle(p->pID); int ep_plus = (plus->pID == p->ID)? 1 : -1; Particle *minus = m_ParticleGrid->GetParticle(p->mID); int ep_minus = (minus->pID == p->ID)? 1 : -1; prop_p.R = (plus->R + plus->N * (plus->len * ep_plus) + minus->R + minus->N * (minus->len * ep_minus)); prop_p.R *= 0.5; prop_p.N = plus->R - minus->R; prop_p.N.normalize(); } else if (p->pID != -1) { Particle *plus = m_ParticleGrid->GetParticle(p->pID); int ep_plus = (plus->pID == p->ID)? 1 : -1; prop_p.R = plus->R + plus->N * (plus->len * ep_plus * 2); prop_p.N = plus->N; } else if (p->mID != -1) { Particle *minus = m_ParticleGrid->GetParticle(p->mID); int ep_minus = (minus->pID == p->ID)? 1 : -1; prop_p.R = minus->R + minus->N * (minus->len * ep_minus * 2); prop_p.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 ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->len,p) + - enc->computeExternalEnergy(p->R,p->N,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); + float prob = exp(ex_energy/m_ExTemp+in_energy/m_InTemp)*m_ShiftProb*p_rev/(m_OptShiftProb+m_ShiftProb*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->GetParticle(pnum); EndPoint P; P.p = p; P.ep = (mtrand->frand() > 0.5)? 1 : -1; RemoveAndSaveTrack(P); if (TrackBackup.m_Probability != 0) { MakeTrackProposal(P); - float prob = (TrackProposal.m_Energy-TrackBackup.m_Energy)/T_in ; + float prob = (TrackProposal.m_Energy-TrackBackup.m_Energy)/m_InTemp ; // prob = exp(prob)*(TrackBackup.proposal_probability) // /(TrackProposal.proposal_probability); prob = exp(prob)*(TrackBackup.m_Probability * pow(del_prob,TrackProposal.m_Length)) /(TrackProposal.m_Probability * pow(del_prob,TrackBackup.m_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.m_Length;k++) { m_ParticleGrid->CreateConnection(T.track[k-1].p,T.track[k-1].ep,T.track[k].p,-T.track[k].ep); } } void MetropolisHastingsSampler::RemoveAndSaveTrack(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; TrackBackup.track[cnt] = Current; EndPoint Next; for (;;) { Next.p = 0; if (Current.ep == 1) { if (Current.p->pID != -1) { Next.p = m_ParticleGrid->GetParticle(Current.p->pID); Current.p->pID = -1; m_ParticleGrid->m_NumConnections--; } } else if (Current.ep == -1) { if (Current.p->mID != -1) { Next.p = m_ParticleGrid->GetParticle(Current.p->mID); Current.p->mID = -1; m_ParticleGrid->m_NumConnections--; } } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 3\n"); break; } if (Next.p == 0) // no successor { Next.ep = 0; // mark as empty successor break; } else { if (Next.p->pID == Current.p->ID) { Next.p->pID = -1; Next.ep = 1; } else if (Next.p->mID == Current.p->ID) { Next.p->mID = -1; Next.ep = -1; } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 4\n"); break; } } ComputeEndPointProposalDistribution(Current); AccumProb *= (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.m_Energy = energy; TrackBackup.m_Probability = AccumProb; TrackBackup.m_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.m_Energy = energy; TrackProposal.m_Probability = AccumProb; TrackProposal.m_Length = cnt; // clear labels for (int j = 0; j < TrackProposal.m_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 c102ef9518..409f15ca5a 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h @@ -1,111 +1,108 @@ /*=================================================================== 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" +#include +#include +#include +#include // ITK #include // MISC #include namespace mitk { class MitkDiffusionImaging_EXPORT MetropolisHastingsSampler { public: - typedef itk::Image< float, 3 > ItkFloatImageType; - ParticleGrid* m_ParticleGrid; - const int* datasz; + ParticleGrid* m_ParticleGrid; + const int* datasz; EnergyComputer* enc; - int m_Iterations; - float width; - float height; - float depth; - 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; + int m_Iterations; + float width; + float height; + float depth; + int m_AcceptedProposals; + + float m_InTemp; + float m_ExTemp; + float m_Density; + + float m_BirthProb; + float m_DeathProb; + float m_ShiftProb; + float m_OptShiftProb; + float m_ConnectionProb; 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; + MTRand* mtrand; Track TrackProposal, TrackBackup; SimpSamp simpsamp; MetropolisHastingsSampler(ParticleGrid* grid, MTRand* randGen); 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/mitkParticle.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h index 127a24e9fc..9b3d15a980 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h @@ -1,79 +1,78 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _PARTICLE #define _PARTICLE #include #include namespace mitk { class MitkDiffusionImaging_EXPORT Particle { public: Particle() { label = 0; pID = -1; mID = -1; inserted = false; } ~Particle() { } vnl_vector_fixed R; vnl_vector_fixed N; - float cap; float len; int gridindex; // index in the grid where it is living int ID; int pID; int mID; int label; int numerator; bool inserted; }; class MitkDiffusionImaging_EXPORT EndPoint { public: EndPoint() {} EndPoint(Particle *p,int ep) { this->p = p; this->ep = ep; } Particle *p; int ep; inline bool operator==(EndPoint P) { return (P.p == p) && (P.ep == ep); } }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp index f3c6fa2d6b..ad1514abaf 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp @@ -1,386 +1,385 @@ /*=================================================================== 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 -#include using namespace mitk; ParticleGrid::ParticleGrid(ItkFloatImageType* image, float cellSize) { // initialize counters m_NumParticles = 0; m_NumConnections = 0; m_NumCellOverflows = 0; // define isotropic grid from voxel spacing and particle length (cellSize) m_GridSize[0] = image->GetLargestPossibleRegion().GetSize()[0]*image->GetSpacing()[0]/cellSize +1; m_GridSize[1] = image->GetLargestPossibleRegion().GetSize()[1]*image->GetSpacing()[1]/cellSize +1; m_GridSize[2] = image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]/cellSize +1; m_GridScale[0] = 1/cellSize; m_GridScale[1] = 1/cellSize; m_GridScale[2] = 1/cellSize; m_CellCapacity = 1024; // maximum number of particles per grid cell m_ContainerCapacity = 100000; // initial particle container capacity int numCells = m_GridSize[0]*m_GridSize[1]*m_GridSize[2]; // number of grid cells m_Particles.resize(m_ContainerCapacity); // allocate and initialize particles m_Grid.resize(numCells*m_CellCapacity, NULL); // allocate and initialize particle grid m_OccupationCount.resize(numCells, 0); // allocate and initialize occupation counter array m_NeighbourTracker.cellidx.resize(8, 0); // allocate and initialize neighbour tracker m_NeighbourTracker.cellidx_c.resize(8, 0); for (int i = 0;i < m_ContainerCapacity;i++) // initialize particle IDs m_Particles[i].ID = i; std::cout << "ParticleGrid: allocated " << (sizeof(Particle)*m_ContainerCapacity + sizeof(Particle*)*m_GridSize[0]*m_GridSize[1]*m_GridSize[2])/1048576 << "mb for " << m_ContainerCapacity/1000 << "k particles." << std::endl; } ParticleGrid::~ParticleGrid() { } bool ParticleGrid::ReallocateGrid() { std::cout << "ParticleGrid: reallocating ..." << std::endl; int new_capacity = m_ContainerCapacity + 100000; // increase container capacity by 100k particles try { m_Particles.resize(new_capacity); // reallocate particles for (int i = 0; i R) { if (m_NumParticles >= m_ContainerCapacity) { if (!ReallocateGrid()) return NULL; } int xint = int(R[0]*m_GridScale[0]); if (xint < 0) return NULL; if (xint >= m_GridSize[0]) return NULL; int yint = int(R[1]*m_GridScale[1]); if (yint < 0) return NULL; if (yint >= m_GridSize[1]) return NULL; int zint = int(R[2]*m_GridScale[2]); if (zint < 0) return NULL; if (zint >= m_GridSize[2]) return NULL; int idx = xint + m_GridSize[0]*(yint + m_GridSize[1]*zint); if (m_OccupationCount[idx] < m_CellCapacity) { Particle *p = &(m_Particles[m_NumParticles]); p->R = R; p->mID = -1; p->pID = -1; m_NumParticles++; p->gridindex = m_CellCapacity*idx + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return p; } else { m_NumCellOverflows++; return NULL; } } bool ParticleGrid::TryUpdateGrid(int k) { Particle* p = &(m_Particles[k]); int xint = int(p->R[0]*m_GridScale[0]); if (xint < 0) return false; if (xint >= m_GridSize[0]) return false; int yint = int(p->R[1]*m_GridScale[1]); if (yint < 0) return false; if (yint >= m_GridSize[1]) return false; int zint = int(p->R[2]*m_GridScale[2]); if (zint < 0) return false; if (zint >= m_GridSize[2]) return false; int idx = xint + m_GridSize[0]*(yint+ zint*m_GridSize[1]); int cellidx = p->gridindex/m_CellCapacity; if (idx != cellidx) // cell has changed { if (m_OccupationCount[idx] < m_CellCapacity) { // remove from old position in grid; int grdindex = p->gridindex; m_Grid[grdindex] = m_Grid[cellidx*m_CellCapacity + m_OccupationCount[cellidx]-1]; m_Grid[grdindex]->gridindex = grdindex; m_OccupationCount[cellidx]--; // insert at new position in grid p->gridindex = idx*m_CellCapacity + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return true; } else { m_NumCellOverflows++; return false; } } return true; } void ParticleGrid::RemoveParticle(int k) { Particle* p = &(m_Particles[k]); int gridIndex = p->gridindex; int cellIdx = gridIndex/m_CellCapacity; int idx = gridIndex%m_CellCapacity; // remove pending connections if (p->mID != -1) DestroyConnection(p,-1); if (p->pID != -1) DestroyConnection(p,+1); // remove from grid if (idx < m_OccupationCount[cellIdx]-1) { m_Grid[gridIndex] = m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1]; m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1] = NULL; m_Grid[gridIndex]->gridindex = gridIndex; } m_OccupationCount[cellIdx]--; // remove from container if (k < m_NumParticles-1) { Particle* last = &m_Particles[m_NumParticles-1]; // last particle // update connections of last particle because its index is changing if (last->mID!=-1) { if ( m_Particles[last->mID].mID == m_NumParticles-1 ) m_Particles[last->mID].mID = k; else if ( m_Particles[last->mID].pID == m_NumParticles-1 ) m_Particles[last->mID].pID = k; } if (last->pID!=-1) { if ( m_Particles[last->pID].mID == m_NumParticles-1 ) m_Particles[last->pID].mID = k; else if ( m_Particles[last->pID].pID == m_NumParticles-1 ) m_Particles[last->pID].pID = k; } m_Particles[k] = m_Particles[m_NumParticles-1]; // move very last particle to empty slot m_Particles[m_NumParticles-1].ID = m_NumParticles-1; // update ID of removed particle to match the index m_Particles[k].ID = k; // update ID of moved particle m_Grid[m_Particles[k].gridindex] = &m_Particles[k]; // update address of moved particle } m_NumParticles--; } void ParticleGrid::ComputeNeighbors(vnl_vector_fixed &R) { float xfrac = R[0]*m_GridScale[0]; float yfrac = R[1]*m_GridScale[1]; float zfrac = R[2]*m_GridScale[2]; 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 >= m_GridSize[0]-1) { xint = m_GridSize[0]-1; dx = -1; } int dy = -1; if (yfrac-yint > 0.5) dy = 1; if (yint <= 0) {yint = 0; dy = 1; } if (yint >= m_GridSize[1]-1) {yint = m_GridSize[1]-1; dy = -1;} int dz = -1; if (zfrac-zint > 0.5) dz = 1; if (zint <= 0) {zint = 0; dz = 1; } if (zint >= m_GridSize[2]-1) {zint = m_GridSize[2]-1; dz = -1;} m_NeighbourTracker.cellidx[0] = xint + m_GridSize[0]*(yint+zint*m_GridSize[1]); m_NeighbourTracker.cellidx[1] = m_NeighbourTracker.cellidx[0] + dx; m_NeighbourTracker.cellidx[2] = m_NeighbourTracker.cellidx[1] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[3] = m_NeighbourTracker.cellidx[2] - dx; m_NeighbourTracker.cellidx[4] = m_NeighbourTracker.cellidx[0] + dz*m_GridSize[0]*m_GridSize[1]; m_NeighbourTracker.cellidx[5] = m_NeighbourTracker.cellidx[4] + dx; m_NeighbourTracker.cellidx[6] = m_NeighbourTracker.cellidx[5] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[7] = m_NeighbourTracker.cellidx[6] - dx; m_NeighbourTracker.cellidx_c[0] = m_CellCapacity*m_NeighbourTracker.cellidx[0]; m_NeighbourTracker.cellidx_c[1] = m_CellCapacity*m_NeighbourTracker.cellidx[1]; m_NeighbourTracker.cellidx_c[2] = m_CellCapacity*m_NeighbourTracker.cellidx[2]; m_NeighbourTracker.cellidx_c[3] = m_CellCapacity*m_NeighbourTracker.cellidx[3]; m_NeighbourTracker.cellidx_c[4] = m_CellCapacity*m_NeighbourTracker.cellidx[4]; m_NeighbourTracker.cellidx_c[5] = m_CellCapacity*m_NeighbourTracker.cellidx[5]; m_NeighbourTracker.cellidx_c[6] = m_CellCapacity*m_NeighbourTracker.cellidx[6]; m_NeighbourTracker.cellidx_c[7] = m_CellCapacity*m_NeighbourTracker.cellidx[7]; m_NeighbourTracker.cellcnt = 0; m_NeighbourTracker.pcnt = 0; } Particle* ParticleGrid::GetNextNeighbor() { if (m_NeighbourTracker.pcnt < m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]]) { return m_Grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt] + (m_NeighbourTracker.pcnt++)]; } else { for(;;) { m_NeighbourTracker.cellcnt++; if (m_NeighbourTracker.cellcnt >= 8) return 0; if (m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]] > 0) break; } m_NeighbourTracker.pcnt = 1; return m_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_Particles[P1->pID]; P1->pID = -1; } else { P2 = &m_Particles[P1->mID]; P1->mID = -1; } if (P2->mID == P1->ID) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } bool ParticleGrid::CheckConsistency() { for (int i=0; iID != i) { std::cout << "Particle ID error!" << std::endl; return false; } if (p->mID!=-1) { Particle* p2 = &m_Particles[p->mID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } if (p->pID!=-1) { Particle* p2 = &m_Particles[p->pID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } } return true; } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h index 724df745f0..281e526bb1 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h @@ -1,119 +1,119 @@ /*=================================================================== 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" +#include // ITK #include namespace mitk { class MitkDiffusionImaging_EXPORT ParticleGrid { public: typedef itk::Image< float, 3 > ItkFloatImageType; int m_NumParticles; // number of particles int m_NumConnections; // number of connections int m_NumCellOverflows; // number of cell overflows ParticleGrid(ItkFloatImageType* image, float cellSize); ~ParticleGrid(); Particle* GetParticle(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); bool CheckConsistency(); protected: bool ReallocateGrid(); std::vector< Particle* > m_Grid; // the grid std::vector< Particle > m_Particles; // particle container std::vector< int > m_OccupationCount; // number of particles per grid cell int m_ContainerCapacity; // maximal number of particles vnl_vector_fixed< int, 3 > m_GridSize; // grid dimensions vnl_vector_fixed< float, 3 > m_GridScale; // scaling factor for grid int m_CellCapacity; // particle capacity of single cell in grid struct NeighborTracker // to run over the neighbors { std::vector< int > cellidx; std::vector< int > cellidx_c; int cellcnt; int pcnt; } m_NeighbourTracker; }; class MitkDiffusionImaging_EXPORT Track { public: std::vector< EndPoint > track; float m_Energy; float m_Probability; int m_Length; Track() { track.resize(1000); } ~Track(){} void clear() { m_Length = 0; m_Energy = 0; m_Probability = 1; } bool isequal(Track& t) { for (int i = 0; i < m_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/GibbsTracking/mitkSimpSamp.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSimpSamp.h index 32e8a31854..acf68a767d 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSimpSamp.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSimpSamp.h @@ -1,121 +1,121 @@ /*=================================================================== 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 _SIMPSAMP #define _SIMPSAMP #include +#include #include -#include "mitkParticle.h" using namespace std; namespace mitk { class MitkDiffusionImaging_EXPORT SimpSamp { float *P; int cnt; public: EndPoint* objs; SimpSamp() { P = (float*) malloc(sizeof(float)*1000); objs = (EndPoint*) malloc(sizeof(EndPoint)*1000); } ~SimpSamp() { free(P); free(objs); } inline void clear() { cnt = 1; P[0] = 0; } inline void add(float p, EndPoint obj) { P[cnt] = P[cnt-1] + p; objs[cnt-1] = obj; cnt++; } inline int draw(float prob) { float r = prob*P[cnt-1]; int j; int rl = 1; int rh = cnt-1; while(rh != rl) { j = rl + (rh-rl)/2; if (r < P[j]) { rh = j; continue; } if (r > P[j]) { rl = j+1; continue; } break; } return rh-1; } inline EndPoint drawObj(float prob) { return objs[draw(prob)]; } inline bool isempty() { if (cnt == 1) return true; else return false; } float probFor(int idx) { return (P[idx+1]-P[idx])/P[cnt-1]; } float probFor(EndPoint& t) { for (int i = 1; i< cnt;i++) { if (t == objs[i-1]) return probFor(i-1); } return 0; } }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSphereInterpolator.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSphereInterpolator.h index f8fc59d1a0..1244715dd8 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSphereInterpolator.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSphereInterpolator.h @@ -1,137 +1,136 @@ /*=================================================================== 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 _SPHEREINTERPOLATOR #define _SPHEREINTERPOLATOR #include class MitkDiffusionImaging_EXPORT SphereInterpolator { public: float *barycoords; int *indices; int size; // size of LUT int sN; // (sizeofLUT-1)/2 int nverts; // number of data vertices float beta; float inva; float b; int *idx; float *interpw; SphereInterpolator(float *barycoords, int *indices, int numverts, int sizeLUT, float beta) { this->barycoords = barycoords; this->indices = indices; this->size = sizeLUT; this->sN = (sizeLUT-1)/2; this->nverts = numverts; this->beta = beta; inva = (sqrt(1+beta)-sqrt(beta)); b = 1/(1-sqrt(1/beta + 1)); } - inline void getInterpolation(vnl_vector_fixed N) { float nx = N[0]; float ny = N[1]; float nz = N[2]; if (nz > 0.5) { int x = float2int(nx); int y = float2int(ny); int i = 3*6*(x+y*size); // (:,1,x,y) idx = indices+i; interpw = barycoords +i; return; } if (nz < -0.5) { int x = float2int(nx); int y = float2int(ny); int i = 3*(1+6*(x+y*size)); // (:,2,x,y) idx = indices+i; interpw = barycoords +i; return; } if (nx > 0.5) { int z = float2int(nz); int y = float2int(ny); int i = 3*(2+6*(z+y*size)); // (:,2,x,y) idx = indices+i; interpw = barycoords +i; return; } if (nx < -0.5) { int z = float2int(nz); int y = float2int(ny); int i = 3*(3+6*(z+y*size)); // (:,2,x,y) idx = indices+i; interpw = barycoords +i; return; } if (ny > 0) { int x = float2int(nx); int z = float2int(nz); int i = 3*(4+6*(x+z*size)); // (:,1,x,y) idx = indices+i; interpw = barycoords +i; return; } else { int x = float2int(nx); int z = float2int(nz); int i = 3*(5+6*(x+z*size)); // (:,1,x,y) idx = indices+i; interpw = barycoords +i; return; } } inline float invrescale(float f) { float x = (fabs(f)-b)*inva; if (f>0) return (x*x-beta); else return beta - x*x; } inline int float2int(float x) { return int((invrescale(x)+1)*sN-0.5); } }; #endif diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index b77fcbc47e..57806429a3 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,389 +1,392 @@ /*=================================================================== 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" +// MITK +#include +#include +#include +#include +#include +#include +#include +#include +// ITK +#include #pragma GCC visibility push(default) #include #pragma GCC visibility pop -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include +// MISC #include +#include +#include namespace itk{ template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::GibbsTrackingFilter(): m_TempStart(0.1), m_TempEnd(0.001), m_NumIt(500000), m_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_BuildFibers(false), m_Steps(10), m_Memory(0), m_ProposalAcceptance(0), m_CurvatureHardThreshold(0.7), m_Meanval_sq(0.0), m_DuplicateImage(true) { } template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::~GibbsTrackingFilter() { } // fill output fiber bundle datastructure template< class ItkQBallImageType > typename GibbsTrackingFilter< ItkQBallImageType >::FiberPolyDataType GibbsTrackingFilter< ItkQBallImageType >::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } return m_FiberPolyData; } template< class ItkQBallImageType > bool GibbsTrackingFilter< ItkQBallImageType > ::EstimateParticleWeight() { MITK_INFO << "itkGibbsTrackingFilter: estimating particle weight"; typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(m_QBallImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); ItkFloatImageType::Pointer gfaImage = gfaFilter->GetOutput(); float samplingStart = 1.0; float samplingStop = 0.66; // GFA iterator typedef ImageRegionIterator< ItkFloatImageType > GfaIteratorType; GfaIteratorType gfaIt(gfaImage, gfaImage->GetLargestPossibleRegion() ); // Mask iterator typedef ImageRegionConstIterator< ItkFloatImageType > MaskIteratorType; MaskIteratorType mit(m_MaskImage, m_MaskImage->GetLargestPossibleRegion() ); // Input iterator typedef ImageRegionConstIterator< ItkQBallImageType > InputIteratorType; InputIteratorType it(m_QBallImage, m_QBallImage->GetLargestPossibleRegion() ); float upper = 0; int count = 0; for(float thr=samplingStart; thr>samplingStop; thr-=0.01) { it.GoToBegin(); mit.GoToBegin(); gfaIt.GoToBegin(); while( !gfaIt.IsAtEnd() ) { if(gfaIt.Get()>thr && mit.Get()>0) { itk::OrientationDistributionFunction odf(it.Get().GetDataPointer()); upper += odf.GetMaxValue()-odf.GetMeanValue(); ++count; } ++it; ++mit; ++gfaIt; } } if (count>0) upper /= count; else return false; m_ParticleWeight = upper/6; return true; } // perform global tracking template< class ItkQBallImageType > void GibbsTrackingFilter< ItkQBallImageType >::GenerateData() { - - m_NumAcceptedFibers = 0; + if (m_QBallImage.IsNull() && m_TensorImage.IsNotNull()) + { + TensorImageToQBallImageFilter::Pointer filter = TensorImageToQBallImageFilter::New(); + filter->SetInput( m_TensorImage ); + filter->Update(); + m_QBallImage = filter->GetOutput(); + } // image sizes and spacing int imgSize[4] = { QBALL_ODFSIZE, m_QBallImage->GetLargestPossibleRegion().GetSize().GetElement(0), m_QBallImage->GetLargestPossibleRegion().GetSize().GetElement(1), m_QBallImage->GetLargestPossibleRegion().GetSize().GetElement(2)}; double spacing[3] = {m_QBallImage->GetSpacing().GetElement(0),m_QBallImage->GetSpacing().GetElement(1),m_QBallImage->GetSpacing().GetElement(2)}; // make sure image has enough slices if (imgSize[1]<3 || imgSize[2]<3 || imgSize[3]<3) { MITK_INFO << "itkGibbsTrackingFilter: image size < 3 not supported"; m_AbortTracking = true; } // calculate rotation matrix vnl_matrix temp = m_QBallImage->GetDirection().GetVnlMatrix(); vnl_matrix directionMatrix; directionMatrix.set_size(3,3); vnl_copy(temp, directionMatrix); vnl_vector_fixed d0 = directionMatrix.get_column(0); d0.normalize(); vnl_vector_fixed d1 = directionMatrix.get_column(1); d1.normalize(); vnl_vector_fixed d2 = directionMatrix.get_column(2); d2.normalize(); directionMatrix.set_column(0, d0); directionMatrix.set_column(1, d1); directionMatrix.set_column(2, d2); vnl_matrix_fixed I = directionMatrix*directionMatrix.transpose(); if(!I.is_identity(mitk::eps)){ MITK_INFO << "itkGibbsTrackingFilter: image direction is not a rotation matrix. Tracking not possible!"; m_AbortTracking = true; } // generate local working copy of QBall image typename ItkQBallImageType::Pointer qballImage; if (m_DuplicateImage) { typedef itk::ImageDuplicator< ItkQBallImageType > DuplicateFilterType; typename DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); duplicator->SetInputImage( m_QBallImage ); duplicator->Update(); qballImage = duplicator->GetOutput(); } else { qballImage = m_QBallImage; } // perform mean subtraction on odfs typedef ImageRegionIterator< ItkQBallImageType > InputIteratorType; InputIteratorType it(qballImage, qballImage->GetLargestPossibleRegion() ); it.GoToBegin(); while (!it.IsAtEnd()) { itk::OrientationDistributionFunction odf(it.Get().GetDataPointer()); float mean = odf.GetMeanValue(); odf -= mean; it.Set(odf.GetDataPointer()); ++it; } // mask image int maskImageSize[3]; float *mask; if(m_MaskImage.IsNotNull()) { mask = (float*) m_MaskImage->GetBufferPointer(); maskImageSize[0] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(0); maskImageSize[1] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(1); maskImageSize[2] = m_MaskImage->GetLargestPossibleRegion().GetSize().GetElement(2); } else { mask = 0; maskImageSize[0] = imgSize[1]; maskImageSize[1] = imgSize[2]; maskImageSize[2] = imgSize[3]; } 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(); - } - else - { - MITK_INFO << "itkGibbsTrackingFilter: unable to open barycoords file"; - m_AbortTracking = true; - } - - ifstream Indices; - lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTIndices.bin"); - Indices.open(lutPath.c_str(), ios::in | ios::binary); - int* ind; - if (Indices.is_open()) - { - int tmp; - ind = new int [1630818]; - Indices.seekg (0, ios::beg); - for (int i=0; i<1630818; i++) - { - Indices.read((char *)&tmp, 4); - ind[i] = tmp; - } - Indices.close(); - } - else - { - MITK_INFO << "itkGibbsTrackingFilter: unable to open indices file"; - m_AbortTracking = true; - } - - // initialize sphere interpolator with lookuptables - SphereInterpolator *sinterp = new SphereInterpolator(coords, ind, QBALL_ODFSIZE, 301, 0.5); - // get paramters float minSpacing; if(spacing[0]m_NumIt) { MITK_INFO << "itkGibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } if (m_CurvatureHardThreshold < mitk::eps) m_CurvatureHardThreshold = 0; unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); - MTRand randGen; -// srand(1); - - MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; + MTRand randGen(1); + srand(1); + SphereInterpolator* interpolator = LoadSphereInterpolator(); + MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; 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); + EnergyComputer encomp(&randGen, qballImage, imgSize,spacing,interpolator,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); sampler->SetEnergyComputer(&encomp); sampler->SetParameters(m_TempStart,singleIts,m_ParticleLength,m_CurvatureHardThreshold,m_ChempotParticle); // main loop + m_NumAcceptedFibers = 0; 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))); 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; } } FiberBuilder fiberBuilder(particleGrid, m_MaskImage); m_FiberPolyData = fiberBuilder.iterate(m_FiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); - delete sinterp; - delete coords; - delete ind; + delete interpolator; delete sampler; delete particleGrid; m_AbortTracking = true; m_BuildFibers = false; MITK_INFO << "itkGibbsTrackingFilter: done generate data"; } + +template< class ItkQBallImageType > +SphereInterpolator* GibbsTrackingFilter< ItkQBallImageType >::LoadSphereInterpolator() +{ + QString applicationDir = QCoreApplication::applicationDirPath(); + applicationDir.append("/"); + mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); + applicationDir.append("../"); + mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); + applicationDir.append("../../"); + mitk::StandardFileLocations::GetInstance()->AddDirectoryForSearch( applicationDir.toStdString().c_str(), false ); + + std::string lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTBaryCoords.bin"); + ifstream BaryCoords; + BaryCoords.open(lutPath.c_str(), ios::in | ios::binary); + float* coords; + if (BaryCoords.is_open()) + { + float tmp; + coords = new float [1630818]; + BaryCoords.seekg (0, ios::beg); + for (int i=0; i<1630818; i++) + { + BaryCoords.read((char *)&tmp, sizeof(tmp)); + coords[i] = tmp; + } + BaryCoords.close(); + } + else + { + MITK_INFO << "itkGibbsTrackingFilter: unable to open barycoords file"; + m_AbortTracking = true; + } + + ifstream Indices; + lutPath = mitk::StandardFileLocations::GetInstance()->FindFile("FiberTrackingLUTIndices.bin"); + Indices.open(lutPath.c_str(), ios::in | ios::binary); + int* ind; + if (Indices.is_open()) + { + int tmp; + ind = new int [1630818]; + Indices.seekg (0, ios::beg); + for (int i=0; i<1630818; i++) + { + Indices.read((char *)&tmp, 4); + ind[i] = tmp; + } + Indices.close(); + } + else + { + MITK_INFO << "itkGibbsTrackingFilter: unable to open indices file"; + m_AbortTracking = true; + } + + // initialize sphere interpolator with lookuptables + return new SphereInterpolator(coords, ind, QBALL_ODFSIZE, 301, 0.5); +} + } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h index 77cb01c7d5..2e97bb4a63 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h @@ -1,154 +1,161 @@ /*=================================================================== 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" +// MITK +#include -#include -#include +// ITK +#include +#include +#include + +// VTK #include #include #include #include #include namespace itk{ template< class ItkQBallImageType > class GibbsTrackingFilter : public ProcessObject { public: typedef GibbsTrackingFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; itkNewMacro(Self) itkTypeMacro( GibbsTrackingFilter, ProcessObject ) + typedef Image< DiffusionTensor3D, 3 > ItkTensorImage; typedef typename ItkQBallImageType::Pointer ItkQBallImageTypePointer; typedef Image< float, 3 > ItkFloatImageType; typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; itkSetMacro( TempStart, float ) itkGetMacro( TempStart, float ) itkSetMacro( TempEnd, float ) itkGetMacro( TempEnd, float ) itkSetMacro( NumIt, unsigned long ) itkGetMacro( NumIt, unsigned long ) 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( CurvatureHardThreshold, float) itkGetMacro( CurvatureHardThreshold, float) itkGetMacro(NumParticles, unsigned long) itkGetMacro(NumConnections, unsigned long) itkGetMacro(NumAcceptedFibers, int) itkGetMacro(ProposalAcceptance, float) itkGetMacro(Steps, unsigned int) // input data itkSetMacro(QBallImage, typename ItkQBallImageType::Pointer) itkSetMacro(MaskImage, ItkFloatImageType::Pointer) + itkSetMacro(TensorImage, ItkTensorImage::Pointer) void GenerateData(); virtual void Update(){ this->GenerateData(); } FiberPolyDataType GetFiberBundle(); protected: GibbsTrackingFilter(); virtual ~GibbsTrackingFilter(); bool EstimateParticleWeight(); + SphereInterpolator* LoadSphereInterpolator(); // Input Images typename ItkQBallImageType::Pointer m_QBallImage; typename ItkFloatImageType::Pointer m_MaskImage; + typename ItkTensorImage::Pointer m_TensorImage; // Tracking parameters float m_TempStart; // Start temperature float m_TempEnd; // End temperature unsigned long m_NumIt; // Total number of iterations unsigned long m_CurrentStep; // current tracking step float m_ParticleWeight; // w (unitless) float m_ParticleWidth; //sigma (mm) float m_ParticleLength; // ell (mm) float m_ChempotConnection; // gross L (chemisches potential) float m_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; int m_NumAcceptedFibers; volatile bool m_BuildFibers; unsigned int m_Steps; float m_Memory; float m_ProposalAcceptance; float m_CurvatureHardThreshold; float m_Meanval_sq; bool m_DuplicateImage; FiberPolyDataType m_FiberPolyData; unsigned long m_NumParticles; unsigned long m_NumConnections; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkGibbsTrackingFilter.cpp" #endif #endif diff --git a/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index df82fc2d51..0536570f09 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,238 +1,239 @@ set(CPP_FILES # DicomImport DicomImport/mitkDicomDiffusionImageReader.cpp DicomImport/mitkGroupDiffusionHeadersFilter.cpp DicomImport/mitkDicomDiffusionImageHeaderReader.cpp DicomImport/mitkGEDicomDiffusionImageHeaderReader.cpp DicomImport/mitkPhilipsDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensDicomDiffusionImageHeaderReader.cpp DicomImport/mitkSiemensMosaicDicomDiffusionImageHeaderReader.cpp # DataStructures IODataStructures/mitkDiffusionImagingObjectFactory.cpp # DataStructures -> DWI IODataStructures/DiffusionWeightedImages/mitkDiffusionImageHeaderInformation.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSource.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageReader.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriter.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageIOFactory.cpp IODataStructures/DiffusionWeightedImages/mitkNrrdDiffusionImageWriterFactory.cpp IODataStructures/DiffusionWeightedImages/mitkDiffusionImageSerializer.cpp # DataStructures -> QBall IODataStructures/QBallImages/mitkQBallImageSource.cpp IODataStructures/QBallImages/mitkNrrdQBallImageReader.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriter.cpp IODataStructures/QBallImages/mitkNrrdQBallImageIOFactory.cpp IODataStructures/QBallImages/mitkNrrdQBallImageWriterFactory.cpp IODataStructures/QBallImages/mitkQBallImage.cpp IODataStructures/QBallImages/mitkQBallImageSerializer.cpp # DataStructures -> Tensor IODataStructures/TensorImages/mitkTensorImageSource.cpp IODataStructures/TensorImages/mitkNrrdTensorImageReader.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriter.cpp IODataStructures/TensorImages/mitkNrrdTensorImageIOFactory.cpp IODataStructures/TensorImages/mitkNrrdTensorImageWriterFactory.cpp IODataStructures/TensorImages/mitkTensorImage.cpp IODataStructures/TensorImages/mitkTensorImageSerializer.cpp # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriter.cpp IODataStructures/FiberBundleX/mitkFiberBundleXReader.cpp IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.cpp IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.cpp IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.cpp # DataStructures -> PlanarFigureComposite IODataStructures/PlanarFigureComposite/mitkPlanarFigureComposite.cpp # DataStructures -> Tbss IODataStructures/TbssImages/mitkTbssImageSource.cpp IODataStructures/TbssImages/mitkTbssRoiImageSource.cpp IODataStructures/TbssImages/mitkNrrdTbssImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssImageIOFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageReader.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageIOFactory.cpp IODataStructures/TbssImages/mitkTbssImage.cpp IODataStructures/TbssImages/mitkTbssRoiImage.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssImageWriterFactory.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriter.cpp IODataStructures/TbssImages/mitkNrrdTbssRoiImageWriterFactory.cpp IODataStructures/TbssImages/mitkTbssImporter.cpp # DataStructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.cpp IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.cpp # Rendering Rendering/vtkMaskedProgrammableGlyphFilter.cpp Rendering/mitkCompositeMapper.cpp Rendering/mitkVectorImageVtkGlyphMapper3D.cpp Rendering/vtkOdfSource.cxx Rendering/vtkThickPlane.cxx Rendering/mitkOdfNormalizationMethodProperty.cpp Rendering/mitkOdfScaleByProperty.cpp Rendering/mitkFiberBundleXMapper2D.cpp Rendering/mitkFiberBundleXMapper3D.cpp Rendering/mitkFiberBundleXThreadMonitorMapper3D.cpp Rendering/mitkTbssImageMapper.cpp Rendering/mitkPlanarCircleMapper3D.cpp Rendering/mitkPlanarPolygonMapper3D.cpp Rendering/mitkConnectomicsNetworkMapper3D.cpp # Interactions Interactions/mitkFiberBundleInteractor.cpp # Algorithms Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.cpp Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.cpp Algorithms/mitkTractAnalyzer.cpp # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.cpp Algorithms/Connectomics/mitkConnectomicsHistogramBase.cpp Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.cpp Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.cpp Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.cpp Algorithms/Connectomics/mitkConnectomicsHistogramCache.cpp Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.cpp Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.cpp # Tractography Tractography/GibbsTracking/mitkParticleGrid.cpp Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp Tractography/GibbsTracking/mitkEnergyComputer.cpp Tractography/GibbsTracking/mitkFiberBuilder.cpp # Function Collection mitkDiffusionFunctionCollection.cpp ) set(H_FILES # function Collection mitkDiffusionFunctionCollection.h # Rendering Rendering/mitkDiffusionImageMapper.h Rendering/mitkTbssImageMapper.h Rendering/mitkOdfVtkMapper2D.h Rendering/mitkFiberBundleXMapper3D.h Rendering/mitkFiberBundleXMapper2D.h Rendering/mitkFiberBundleXThreadMonitorMapper3D.h Rendering/mitkPlanarCircleMapper3D.h Rendering/mitkPlanarPolygonMapper3D.h Rendering/mitkConnectomicsNetworkMapper3D.h # Reconstruction Reconstruction/itkDiffusionQballReconstructionImageFilter.h Reconstruction/mitkTeemDiffusionTensor3DReconstructionImageFilter.h Reconstruction/itkAnalyticalDiffusionQballReconstructionImageFilter.h Reconstruction/itkDiffusionMultiShellQballReconstructionImageFilter.h Reconstruction/itkPointShell.h Reconstruction/itkOrientationDistributionFunction.h Reconstruction/itkDiffusionIntravoxelIncoherentMotionReconstructionImageFilter.h Reconstruction/itkRegularizedIVIMLocalVariationImageFilter.h Reconstruction/itkRegularizedIVIMReconstructionFilter.h Reconstruction/itkRegularizedIVIMReconstructionSingleIteration.h # IO Datastructures IODataStructures/DiffusionWeightedImages/mitkDiffusionImage.h IODataStructures/TbssImages/mitkTbssImporter.h # DataStructures -> FiberBundleX IODataStructures/FiberBundleX/mitkFiberBundleX.h IODataStructures/FiberBundleX/mitkFiberBundleXWriter.h IODataStructures/FiberBundleX/mitkFiberBundleXReader.h IODataStructures/FiberBundleX/mitkFiberBundleXIOFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXWriterFactory.h IODataStructures/FiberBundleX/mitkFiberBundleXSerializer.h IODataStructures/FiberBundleX/mitkFiberBundleXThreadMonitor.h # Datastructures Connectomics IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetwork.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkReader.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkIOFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkSerializer.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriter.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkWriterFactory.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsNetworkDefinitions.h IODataStructures/ConnectomicsNetwork/mitkConnectomicsConstantsManager.h # Tractography Tractography/itkGibbsTrackingFilter.h Tractography/itkStochasticTractographyFilter.h Tractography/itkStreamlineTrackingFilter.h Tractography/GibbsTracking/mitkParticle.h Tractography/GibbsTracking/mitkParticleGrid.h Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h Tractography/GibbsTracking/mitkSimpSamp.h Tractography/GibbsTracking/mitkEnergyComputer.h Tractography/GibbsTracking/mitkSphereInterpolator.h Tractography/GibbsTracking/mitkFiberBuilder.h # Algorithms Algorithms/itkDiffusionQballGeneralizedFaImageFilter.h Algorithms/itkDiffusionQballPrepareVisualizationImageFilter.h Algorithms/itkTensorDerivedMeasurementsFilter.h Algorithms/itkBrainMaskExtractionImageFilter.h Algorithms/itkB0ImageExtractionImageFilter.h Algorithms/itkB0ImageExtractionToSeparateImageFilter.h Algorithms/itkTensorImageToDiffusionImageFilter.h Algorithms/itkTensorToL2NormImageFilter.h Algorithms/itkTractDensityImageFilter.h Algorithms/itkTractsToFiberEndingsImageFilter.h Algorithms/itkTractsToRgbaImageFilter.h Algorithms/itkGaussianInterpolateImageFunction.h Algorithms/mitkPartialVolumeAnalysisHistogramCalculator.h Algorithms/mitkPartialVolumeAnalysisClusteringCalculator.h Algorithms/itkDiffusionTensorPrincipleDirectionImageFilter.h Algorithms/itkCartesianToPolarVectorImageFilter.h Algorithms/itkPolarToCartesianVectorImageFilter.h Algorithms/itkDistanceMapFilter.h Algorithms/itkProjectionFilter.h Algorithms/itkSkeletonizationFilter.h Algorithms/itkReduceDirectionGradientsFilter.h Algorithms/itkResidualImageFilter.h Algorithms/itkExtractChannelFromRgbaImageFilter.h + Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h # Algorithms Connectomics Algorithms/Connectomics/mitkConnectomicsNetworkCreator.h Algorithms/Connectomics/mitkConnectomicsHistogramBase.h Algorithms/Connectomics/mitkConnectomicsDegreeHistogram.h Algorithms/Connectomics/mitkConnectomicsShortestPathHistogram.h Algorithms/Connectomics/mitkConnectomicsBetweennessHistogram.h Algorithms/Connectomics/mitkConnectomicsHistogramCache.h Algorithms/Connectomics/mitkConnectomicsSyntheticNetworkGenerator.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingPermutationModularity.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingManager.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionBase.h Algorithms/Connectomics/mitkConnectomicsSimulatedAnnealingCostFunctionModularity.h - Algorithms/itkTensorReconstructionWithEigenvalueCorrectionFilter.h + ) set( TOOL_FILES ) if(WIN32) endif(WIN32) #MITK_MULTIPLEX_PICTYPE( Algorithms/mitkImageRegistrationMethod-TYPE.cpp )