diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp deleted file mode 100644 index 442c343801..0000000000 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/BuildFibres.cpp +++ /dev/null @@ -1,194 +0,0 @@ -/*=================================================================== - -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 "matrix.h" -#define _USE_MATH_DEFINES -#include - -#include -#include -#include - -using namespace std; - -#define PI M_PI - -#include "mitkParticleGrid.h" - -#include -#include -#include -#include -#include -#include - -#include -#include - -class FiberBuilder -{ -public: - Particle *particles; - int pcnt; - int attrcnt; - typedef vector< Particle* > ParticleContainerType; - typedef vector< ParticleContainerType* > FiberContainerType; - - vtkSmartPointer m_VtkCellArray; - vtkSmartPointer m_VtkPoints; - - typedef itk::Vector OdfVectorType; - typedef itk::Image ItkQBallImgType; - ItkQBallImgType::Pointer m_ItkQBallImage; - float m_FiberLength; - itk::Point m_LastPoint; - - FiberBuilder(float *points, int numPoints, double spacing[], ItkQBallImgType::Pointer 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 = pVector(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 = pVector(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() - { - free(particles); - } - - vtkSmartPointer iterate(int minFiberLength) - { - int cur_label = 1; - int numFibers = 0; - m_FiberLength = 0; - for (int k = 0; k < pcnt;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; - } - } - vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); - fiberPolyData->SetPoints(m_VtkPoints); - fiberPolyData->SetLines(m_VtkCellArray); -// vtkSmartPointer cleaner = vtkSmartPointer::New(); -// cleaner->SetInput(fiberPolyData); -// cleaner->Update(); -// fiberPolyData = cleaner->GetOutput(); - return fiberPolyData; - } - - void 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]; - itk::Point point; - m_ItkQBallImage->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 labelPredecessors(Particle *dp, vtkSmartPointer container) - { - if (dp->mID != -1 && dp->mID!=dp->ID) - { - if (dp->ID!=particles[dp->mID].pID) - { - if (dp->ID==particles[dp->mID].mID) - { - int tmp = particles[dp->mID].pID; - particles[dp->mID].pID = particles[dp->mID].mID; - particles[dp->mID].mID = tmp; - } - } - if (particles[dp->mID].label == 0) - { - particles[dp->mID].label = dp->label; - particles[dp->mID].numerator = dp->numerator-1; - labelPredecessors(&(particles[dp->mID]), container); - } - } - - AddPoint(dp, container); - } - - void 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[dp->pID].pID) - { - int tmp = particles[dp->pID].pID; - particles[dp->pID].pID = particles[dp->pID].mID; - particles[dp->pID].mID = tmp; - } - } - if (particles[dp->pID].label == 0) - { - particles[dp->pID].label = dp->label; - particles[dp->pID].numerator = dp->numerator+1; - labelSuccessors(&(particles[dp->pID]), container); - } - } - } -}; - -#endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/EnergyComputer_connec.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/EnergyComputer_connec.cpp deleted file mode 100644 index 339e5d859c..0000000000 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/EnergyComputer_connec.cpp +++ /dev/null @@ -1,429 +0,0 @@ -/*=================================================================== - -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 "mitkParticleGrid.h" -#include "SphereInterpolator.cpp" -#include -#include - -using namespace mitk; - -inline float 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; -} - -inline float 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); - -} - -class EnergyComputer -{ - -public: - - - float eigen_energy; - vnl_matrix_fixed m_RotationMatrix; - float *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; - - - 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(float *data, const int *dsz, double *cellsize, SphereInterpolator *sp, ParticleGrid *pcon, float *spimg, int spmult, vnl_matrix_fixed rotMatrix) - { - 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 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 drawSpatPosition(pVector *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->SetXYZ(voxsize_sp_w*((float)(spatidx[rh-1] % w_sp) + mtrand.frand()), - voxsize_sp_h*((float)((spatidx[rh-1]/w_sp) % h_sp) + mtrand.frand()), - voxsize_sp_d*((float)(spatidx[rh-1]/(w_sp*h_sp)) + mtrand.frand())); - } - - float SpatProb(pVector R) - { - int rx = int(R.GetX()/voxsize_sp_w); - int ry = int(R.GetY()/voxsize_sp_h); - int rz = int(R.GetZ()/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; - } - - inline float evaluateODF(pVector &R, pVector &N, float &len) - { - const int CU = 10; - pVector 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; - - 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; - - 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; - - 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; - - 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; - - 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; - - 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; - - 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; - - 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; - - - } - - - } - - Dn *= 1.0/(2*CU+1); - return Dn; - } - - - //////////////////////////////////////////////////////////////////////////// - ////// External Energy - //////////////////////////////////////////////////////////////////////////// - inline float computeExternalEnergy(pVector &R, pVector &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(N*p->N); - float bw = mbesseli0(dot); - float dpos = (p->R-R).norm_square(); - float w = mexp(dpos*gamma_s); - Sn += w*(bw+chempot2)*p->cap ; - //Sn += w*(bw-meanval_sq)*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; - //energy += (2*(Dn/particle_weight-Sn) - (mbesseli0(1.0)-meanval_sq)*cap)*cap; - - return energy*ex_strength; - } - - - //////////////////////////////////////////////////////////////////////////// - ////// Internal Energy - //////////////////////////////////////////////////////////////////////////// - - inline float computeInternalEnergy(Particle *dp) - { - float energy = eigen_energy; - - if (dp->pID != -1) - energy += computeInternalEnergyConnection(dp,+1); - - if (dp->mID != -1) - energy += computeInternalEnergyConnection(dp,-1); - - //ie_file << energy << "\n"; - - return energy; - } - - inline float 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); - } - - inline float computeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) - { -#ifdef TIMING - tic(&internalenergy_time); -#endif - - if ((p1->N*p2->N)*ep1*ep2 > -curv_hard) - return -INFINITY; - - pVector R1 = p1->R + (p1->N * (p1->len * ep1)); - pVector R2 = p2->R + (p2->N * (p2->len * ep2)); - - if ((R1-R2).norm_square() > particle_length_sq) - return -INFINITY; - - pVector R = (p2->R + p1->R)*0.5; - - if (SpatProb(R) == 0) - return -INFINITY; - - float norm1 = (R1-R).norm_square(); - float norm2 = (R2-R).norm_square(); - - - float energy = (eigencon_energy-norm1-norm2)*in_strength; - -#ifdef TIMING - toc(&internalenergy_time); -#endif - - return energy; - } -}; - -#endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMCBase.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMCBase.cpp deleted file mode 100644 index 48b34fe294..0000000000 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMCBase.cpp +++ /dev/null @@ -1,153 +0,0 @@ -/*=================================================================== - -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 "EnergyComputer_connec.cpp" -#include - -class RJMCMCBase -{ -public: - - 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; - - RJMCMCBase(float *points,int numPoints, float *dimg, const int *dsz, double *voxsize, double cellsize) - : m_QBallImageData(dimg) - , datasz(dsz) - , enc(0) - , width(dsz[1]*voxsize[0]) - , height(dsz[2]*voxsize[1]) - , depth(dsz[3]*voxsize[2]) - , voxsize(voxsize) - , m_NumAttributes(0) - , m_AcceptedProposals(0) - { - m_ParticleGrid = new ParticleGrid(); - 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,"RJMCMCBase: out of Memory!\n"); - return; - } - - m_NumAttributes = 10; - for (int k = 0; k < numPoints; k++) - { - Particle *p = m_ParticleGrid->newParticle(pVector(points[m_NumAttributes*k], points[m_NumAttributes*k+1],points[m_NumAttributes*k+2])); - if (p!=0) - { - p->N = pVector(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_Iterations = 0; - m_AcceptedProposals = 0; - } - - ~RJMCMCBase() - { - delete m_ParticleGrid; - } - - void 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.GetX(); - npoints[m_NumAttributes*k+1] = p->R.GetY(); - npoints[m_NumAttributes*k+2] = p->R.GetZ(); - npoints[m_NumAttributes*k+3] = p->N.GetX(); - npoints[m_NumAttributes*k+4] = p->N.GetY(); - npoints[m_NumAttributes*k+5] = p->N.GetZ(); - 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); - } - } - - void SetEnergyComputer(EnergyComputer *e) - { - enc = e; - } - - void 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; - } - - virtual void IterateOneStep() - { - - } -}; - - - diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMC_randshift.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMC_randshift.cpp deleted file mode 100644 index c5172e5ac1..0000000000 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/RJMCMC_randshift.cpp +++ /dev/null @@ -1,625 +0,0 @@ -/*=================================================================== - -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 "RJMCMCBase.cpp" -#include - -class RJMCMC : public RJMCMCBase -{ -public: - - 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; - - Track TrackProposal, TrackBackup; - SimpSamp simpsamp; - - RJMCMC(float *points,int numPoints, float *dimg, const int *dsz, double *voxsz, double cellsz) : RJMCMCBase(points,numPoints,dimg,dsz,voxsz,cellsz) - { - externalEnergy = 0; - internalEnergy = 0; - } - - void 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 SetTemperature(float temp) - { - T_in = temp; - dens = exp(-m_ChempotParticle/T_in); - } - - void IterateOneStep() - { - float randnum = mtrand.frand(); - //randnum = 0; - - /////////////////////////////////////////////////////////////// - //////// Birth Proposal - /////////////////////////////////////////////////////////////// - if (randnum < p_birth) - { - -#ifdef TIMING - tic(&birthproposal_time); - birthstats.propose(); -#endif - - pVector 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); - - pVector N; 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; - - prop_p.R.distortn(sigma_g); - prop_p.N.distortn(sigma_g/(2*p->len)); - 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) - { - pVector Rtmp = p->R; - pVector 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))*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 = prop_p.N*p->N; - float p_rev = exp(-((prop_p.R-p->R).norm_square() + (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) - { - pVector Rtmp = p->R; - pVector 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 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 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,"RJMCMC_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,"RJMCMC_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 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(); - - // 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 ComputeEndPointProposalDistribution(EndPoint P) - { - Particle *p = P.p; - int ep = P.ep; - - float dist,dot; - pVector 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).norm_square(); - if (dist < dthres) - { - dot = 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).norm_square(); - if (dist < dthres) - { - dot = 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/auxilary_classes.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/auxilary_classes.cpp deleted file mode 100644 index e1f6dedd8f..0000000000 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/auxilary_classes.cpp +++ /dev/null @@ -1,514 +0,0 @@ -/*=================================================================== - -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 _AUXCLASS -#define _AUXCLASS - -#ifndef INFINITY -#define INFINITY 1000000000 -#endif - -#include "MersenneTwister.h" -#include -using namespace std; - -MTRand mtrand; - -class pVector -{ -private: - float x; - float y; - float z; - -public: - - pVector() - { - - } - - pVector(float x,float y,float z) - { - this->x = x; - this->y = y; - this->z = z; - } - - inline void SetXYZ(float sx,float sy, float sz) - { - x = sx; - y = sy; - z = sz; - } - - inline void GetXYZ(float *xyz) - { - xyz[0] = x; - xyz[1] = y; - xyz[2] = z; - } - - inline float GetX() - { - return x; - } - - inline float GetY() - { - return y; - } - - inline float GetZ() - { - return z; - } - - inline void rand(float w, float h, float d) - { - this->x = mtrand.frand()*w; - this->y = mtrand.frand()*h; - this->z = mtrand.frand()*d; - } - - inline void rand_sphere() - { - this->x = mtrand.frandn(); - this->y = mtrand.frandn(); - this->z = mtrand.frandn(); - normalize(); - } - - inline void normalize() - { - float norm = sqrt(x*x+y*y+z*z)+ 0.00000001; - *this /= norm; - } - - inline float norm_square() - { - return x*x + y*y + z*z; - } - - inline void distortn(float sigma) - { - x += sigma*mtrand.frandn(); - y += sigma*mtrand.frandn(); - z += sigma*mtrand.frandn(); - } - - inline float operator[](int index) - { - switch(index) - { - case 0: - return x; - case 1: - return y; - case 2: - return z; - default: - return 0.0f; - } - } - - inline pVector operator*(float s) - { - return pVector(s*x,s*y,s*z); - } - - inline void operator*=(float &s) - { - x *= s; - y *= s; - z *= s; - } - - inline pVector operator+(pVector R) - { - return pVector(x+R.x,y+R.y,z+R.z); - } - - inline void operator+=(pVector R) - { - x += R.x; - y += R.y; - z += R.z; - } - - inline pVector operator-(pVector R) - { - return pVector(x-R.x,y-R.y,z-R.z); - } - - inline void operator-=(pVector R) - { - x -= R.x; - y -= R.y; - z -= R.z; - } - - inline pVector operator/(float &s) - { - return pVector(x/s,y/s,z/s); - } - - inline void operator/=(float &s) - { - x /= s; - y /= s; - z /= s; - } - - - inline float operator*(pVector R) - { - return x*R.x+y*R.y+z*R.z; - } -}; - -class Particle -{ -public: - - Particle() - { - label = 0; - pID = -1; - mID = -1; - inserted = false; - } - - ~Particle() - { - } - - pVector R; - pVector 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 EnergyGradient -{ -public: - pVector gR; - pVector gN; - - inline float norm2() - { - return gR.norm_square() + gN.norm_square(); - } -} ; - - -template -class SimpSamp -{ - - float *P; - int cnt; - -public: - T *objs; - - - SimpSamp() - { - P = (float*) malloc(sizeof(float)*1000); - objs = (T*) malloc(sizeof(T)*1000); - } - ~SimpSamp() - { - free(P); - free(objs); - } - - inline void clear() - { - cnt = 1; - P[0] = 0; - } - - inline void add(float p, T obj) - { - P[cnt] = P[cnt-1] + p; - objs[cnt-1] = obj; - cnt++; - } - - // inline int draw() - // { - // float r = mtrand.frand()*P[cnt-1]; - // for (int i = 1; i < cnt; i++) - // { - // if (r <= P[i]) - // return i-1; - // } - // return cnt-2; - // } - - inline int draw() - { - float r = mtrand.frand()*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 T drawObj() - { - return objs[draw()]; - } - - 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(T &t) - { - for (int i = 1; i< cnt;i++) - { - if (t == objs[i-1]) - return probFor(i-1); - } - return 0; - } - - - -}; - - -class 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); - } -}; - -class 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; - } - -}; - -float getMax(float *arr, int cnt) -{ - float max = arr[0]; - for (int i = 1; i < cnt; i++) - { - if (arr[i] > max) - max = arr[i]; - } - return max; -} - - - -float getMin(float *arr, int cnt) -{ - float min = arr[0]; - for (int i = 1; i < cnt; i++) - { - if (arr[i] < min) - min = arr[i]; - } - return min; -} - - -int getArgMin(float *arr, int cnt) -{ - float min = arr[0]; - int idx = 0; - for (int i = 1; i < cnt; i++) - { - if (arr[i] < min) - { - min = arr[i]; - idx = i; - } - } - return idx; -} - - - - -inline float distLseg(pVector &R1,pVector &N1,pVector &R2,pVector &N2,float &len) -{ - - pVector D = R1-R2; - float beta = N1*N2; - float divisor = 1.001-beta*beta; - float gamma1 = N1*D; - float gamma2 = N2*D; - float t,u; - float EPdist[4]; - - pVector Q; - float dist = 102400000.0; - - while(true) - { - - t = -(gamma1+beta*gamma2) / divisor; - u = (gamma1*beta+gamma2) / divisor; - if (fabs(t) < len && fabs(u) < len) - { - Q = D + N1*t - N2*u; - dist = Q*Q; - break; - } - - beta = len*beta; - - t = beta - gamma1; - if (fabs(t) < len) - { - Q = D + N1*t - N2*len; - float d = Q*Q; - if (d < dist) dist = d; - } - - t = -beta - gamma1; - if (fabs(t) < len) - { - Q = D + N1*t + N2*len; - float d = Q*Q; - if (d < dist) dist = d; - } - - u = beta + gamma2; - if (fabs(u) < len) - { - Q = D + N1*len - N2*u; - float d = Q*Q; - if (d < dist) dist = d; - } - - u = -beta + gamma2; - if (fabs(u) < len) - { - Q = D - N1*len - N2*u; - float d = Q*Q; - if (d < dist) dist = d; - } - - if (dist != 102400000.0) - break; - - - EPdist[0] = beta + gamma1 - gamma2; - EPdist[1] = -beta + gamma1 + gamma2; - EPdist[2] = -beta - gamma1 - gamma2; - EPdist[3] = beta - gamma1 + gamma2; - int c = getArgMin(EPdist,4); - if (c==0) {t = +len; u = +len; } - if (c==1) {t = +len; u = -len; } - if (c==2) {t = -len; u = +len; } - if (c==3) {t = -len; u = -len; } - Q = D + N1*t - N2*u; - dist = Q*Q; - break; - - } - - - return dist; - -} - - -#endif - - diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp new file mode 100644 index 0000000000..ec700a555a --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.cpp @@ -0,0 +1,339 @@ +/*=================================================================== + +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) +{ + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + + 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; + } + } + 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); +} diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h new file mode 100644 index 0000000000..6a90719beb --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkEnergyComputer.h @@ -0,0 +1,97 @@ +/*=================================================================== + +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 "mitkParticleGrid.h" +#include "mitkSphereInterpolator.h" +#include "MersenneTwister.h" +#include +#include + +using namespace mitk; + +class MitkDiffusionImaging_EXPORT EnergyComputer +{ + +public: + float eigen_energy; + vnl_matrix_fixed m_RotationMatrix; + float *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); + + 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 new file mode 100644 index 0000000000..095fa63b58 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp @@ -0,0 +1,147 @@ +/*=================================================================== + +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) +{ + 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) +{ + int cur_label = 1; + int numFibers = 0; + m_FiberLength = 0; + for (int k = 0; k < pcnt;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; + } + } + 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]; + itk::Point point; + m_ItkQBallImage->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[dp->mID].mID) + { + int tmp = particles[dp->mID].pID; + particles[dp->mID].pID = particles[dp->mID].mID; + particles[dp->mID].mID = tmp; + } + } + if (particles[dp->mID].label == 0) + { + particles[dp->mID].label = dp->label; + particles[dp->mID].numerator = dp->numerator-1; + labelPredecessors(&(particles[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[dp->pID].pID) + { + int tmp = particles[dp->pID].pID; + particles[dp->pID].pID = particles[dp->pID].mID; + particles[dp->pID].mID = tmp; + } + } + if (particles[dp->pID].label == 0) + { + particles[dp->pID].label = dp->label; + particles[dp->pID].numerator = dp->numerator+1; + labelSuccessors(&(particles[dp->pID]), container); + } + } +} diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h new file mode 100644 index 0000000000..0f0593095a --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.h @@ -0,0 +1,73 @@ +/*=================================================================== + +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::Vector OdfVectorType; + typedef itk::Image ItkQBallImgType; + ItkQBallImgType::Pointer m_ItkQBallImage; + float m_FiberLength; + itk::Point m_LastPoint; + + FiberBuilder(float *points, int numPoints, double spacing[], ItkQBallImgType::Pointer image); + + ~FiberBuilder(); + + vtkSmartPointer iterate(int minFiberLength); + + void AddPoint(Particle *dp, vtkSmartPointer container); + + void labelPredecessors(Particle *dp, vtkSmartPointer container); + + void labelSuccessors(Particle *dp, vtkSmartPointer container); +}; + +} + +#endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp new file mode 100644 index 0000000000..266a54d546 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp @@ -0,0 +1,704 @@ +/*=================================================================== + +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) + , m_AcceptedProposals(0) +{ + mtrand = rgen; + + 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_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); + } +} + +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 new file mode 100644 index 0000000000..3df94e930a --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h @@ -0,0 +1,107 @@ +/*=================================================================== + +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 + +#include +#include "mitkParticleGrid.h" +#include "mitkEnergyComputer.h" +#include "MersenneTwister.h" +#include "mitkSimpSamp.h" +#include + +namespace mitk +{ + +class MitkDiffusionImaging_EXPORT MetropolisHastingsSampler +{ +public: + + + 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); + + 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/mitkParticle.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h new file mode 100644 index 0000000000..127a24e9fc --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticle.h @@ -0,0 +1,79 @@ +/*=================================================================== + +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 new file mode 100644 index 0000000000..ce7afad1cb --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp @@ -0,0 +1,430 @@ +/*=================================================================== + +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() +{ + //// 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 ParticleGrid::allocate(int _capacity, int _nx, int _ny, int _nz, float cellsize, int cellcapacity) +{ + //// involving the container + capacity = _capacity; + 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 + } + + ////// involvin the grid + nx = _nx; ny = _ny; nz = _nz; csize = cellcapacity; + 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; +} + + + +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) +{ + 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 new file mode 100644 index 0000000000..4c4b48b8bb --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h @@ -0,0 +1,118 @@ +/*=================================================================== + +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 + +#include "mitkParticle.h" +#include "MitkDiffusionImagingExports.h" + +namespace mitk +{ + +class MitkDiffusionImaging_EXPORT ParticleGrid +{ + + //////////////// Container +public: + 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(); + + float GetMemoryUsage(); + int allocate(int _capacity, int _nx, int _ny, int _nz, float cellsize, int cellcapacity); + int reallocate(); + + int ID_2_index(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/GibbsTracking/mitkSimpSamp.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSimpSamp.h new file mode 100644 index 0000000000..32e8a31854 --- /dev/null +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSimpSamp.h @@ -0,0 +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 "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/SphereInterpolator.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSphereInterpolator.h similarity index 86% rename from Modules/DiffusionImaging/Tractography/GibbsTracking/SphereInterpolator.cpp rename to Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSphereInterpolator.h index f13bfbd18b..f8fc59d1a0 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/SphereInterpolator.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkSphereInterpolator.h @@ -1,137 +1,137 @@ /*=================================================================== 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 _SPHEREINTERPOLATOR +#define _SPHEREINTERPOLATOR -#include "auxilary_classes.cpp" +#include - -class SphereInterpolator +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) + 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 3c8588b58c..16c92c3010 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,664 +1,666 @@ /*=================================================================== 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 "itkPointShell.h" -#include "GibbsTracking/BuildFibres.cpp" +#include "GibbsTracking/mitkFiberBuilder.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(): 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) { 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() { 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 > bool GibbsTrackingFilter< TInputOdfImage, TInputROIImage > ::EstimateParticleWeight() { MITK_INFO << "itkGibbsTrackingFilter: estimating particle weight"; typedef itk::DiffusionQballGeneralizedFaImageFilter GfaFilterType; GfaFilterType::Pointer gfaFilter = GfaFilterType::New(); gfaFilter->SetInput(m_ItkQBallImage); gfaFilter->SetComputationMethod(GfaFilterType::GFA_STANDARD); gfaFilter->Update(); m_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< GfaImageType > GfaIteratorType; GfaIteratorType gfaIt(m_GfaImage, m_GfaImage->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() ); float upper = 0; int count = 0; for(float thr=samplingStart; thr>samplingStop; thr-=0.01) { git.GoToBegin(); gfaIt.GoToBegin(); while( !gfaIt.IsAtEnd() ) { if(gfaIt.Get()>thr) { itk::OrientationDistributionFunction odf(git.Get().GetDataPointer()); upper += odf.GetMaxValue()-odf.GetMeanValue(); ++count; } ++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(){ // 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)}; // make sure image has enough slices if (qBallImageSize[1]<3 || qBallImageSize[2]<3 || qBallImageSize[3]<3) { MITK_INFO << "itkGibbsTrackingFilter: image size < 3 not supported"; m_AbortTracking = true; } // calculate rotation matrix - vnl_matrix_fixed directionMatrix = m_ItkQBallImage->GetDirection().GetVnlMatrix(); - 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(); + vnl_matrix temp = m_ItkQBallImage->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(); + 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; i0 && i%qBallImageSize[0] == 0 && i>0) { sum /= qBallImageSize[0]; for (int j=i-qBallImageSize[0]; jGetBufferPointer(); 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]; } int mask_oversamp_mult = maskImageSize[0]/qBallImageSize[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(qBallImageSpacing[0]m_NumIt) { MITK_INFO << "itkGibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } if (m_CurvatureHardThreshold < mitk::eps) m_CurvatureHardThreshold = 0; unsigned long singleIts = (unsigned long)((1.0*m_NumIt) / (1.0*m_Steps)); + MTRand randGen(1); + srand(1); + // setup metropolis hastings sampler MITK_INFO << "itkGibbsTrackingFilter: setting up MH-sampler"; if (m_Sampler!=NULL) delete m_Sampler; - m_Sampler = new RJMCMC(NULL, 0, workingQballImage, qBallImageSize, qBallImageSpacing, cellsize); - - mtrand.seed((unsigned long)0); - srand(0); + m_Sampler = new MetropolisHastingsSampler(&randGen, NULL, 0, workingQballImage, qBallImageSize, qBallImageSpacing, cellsize); // setup energy computer MITK_INFO << "itkGibbsTrackingFilter: setting up Energy-computer"; - EnergyComputer encomp(workingQballImage,qBallImageSize,qBallImageSpacing,sinterp,m_Sampler->m_ParticleGrid,mask,mask_oversamp_mult, directionMatrix); + 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); 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; // 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; } } 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; delete sinterp; delete coords; delete ind; delete workingQballImage; 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 8a5a356806..5919f6e018 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h @@ -1,200 +1,198 @@ /*=================================================================== 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/auxilary_classes.cpp" -#include "GibbsTracking/RJMCMC_randshift.cpp" -#include "GibbsTracking/EnergyComputer_connec.cpp" +#include "GibbsTracking/MersenneTwister.h" +#include "GibbsTracking/mitkMetropolisHastingsSampler.h" +#include "GibbsTracking/mitkEnergyComputer.h" #include #include #include #include #include #include #include -#include - namespace itk{ template< class TInputQBallImage, class TInputROIImage > 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 ); + 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; - itkSetMacro( TempStart, float ); - itkGetMacro( TempStart, float ); + itkSetMacro( TempStart, float ) + itkGetMacro( TempStart, float ) - itkSetMacro( TempEnd, float ); - itkGetMacro( TempEnd, float ); + itkSetMacro( TempEnd, float ) + itkGetMacro( TempEnd, float ) - itkSetMacro( NumIt, unsigned long ); - itkGetMacro( NumIt, unsigned long ); + itkSetMacro( NumIt, unsigned long ) + itkGetMacro( NumIt, unsigned long ) - itkSetMacro( ParticleWeight, float ); - itkGetMacro( ParticleWeight, float ); + itkSetMacro( ParticleWeight, float ) + itkGetMacro( ParticleWeight, float ) /** width of particle sigma (std-dev of gaussian around center) **/ - itkSetMacro( ParticleWidth, float ); - itkGetMacro( ParticleWidth, float ); + itkSetMacro( ParticleWidth, float ) + itkGetMacro( ParticleWidth, float ) /** length of particle from midpoint to ends **/ - itkSetMacro( ParticleLength, float ); - itkGetMacro( ParticleLength, float ); + itkSetMacro( ParticleLength, float ) + itkGetMacro( ParticleLength, float ) - itkSetMacro( ChempotConnection, float ); - itkGetMacro( ChempotConnection, float ); + itkSetMacro( ChempotConnection, float ) + itkGetMacro( ChempotConnection, float ) - itkSetMacro( ChempotParticle, float ); - itkGetMacro( ChempotParticle, float ); + itkSetMacro( ChempotParticle, float ) + itkGetMacro( ChempotParticle, float ) - itkSetMacro( InexBalance, float ); - itkGetMacro( InexBalance, float ); + itkSetMacro( InexBalance, float ) + itkGetMacro( InexBalance, float ) - itkSetMacro( Chempot2, float ); - itkGetMacro( Chempot2, float ); + itkSetMacro( Chempot2, float ) + itkGetMacro( Chempot2, float ) - itkSetMacro( FiberLength, int ); - itkGetMacro( FiberLength, int ); + itkSetMacro( FiberLength, int ) + itkGetMacro( FiberLength, int ) - itkSetMacro( AbortTracking, bool ); - itkGetMacro( AbortTracking, bool ); + itkSetMacro( AbortTracking, bool ) + itkGetMacro( AbortTracking, bool ) - itkSetMacro( CurrentStep, unsigned long ); - itkGetMacro( CurrentStep, unsigned long ); + itkSetMacro( CurrentStep, unsigned long ) + itkGetMacro( CurrentStep, unsigned long ) - itkSetMacro( SubtractMean, bool); - itkGetMacro( SubtractMean, bool); + itkSetMacro( SubtractMean, bool) + itkGetMacro( SubtractMean, bool) - itkSetMacro( CurvatureHardThreshold, float); - itkGetMacro( CurvatureHardThreshold, float); + itkSetMacro( CurvatureHardThreshold, float) + itkGetMacro( CurvatureHardThreshold, float) /** Set/Get the Odf Input Image **/ - itkSetInputMacro(OdfImage, InputQBallImageType, 0); - itkGetInputMacro(OdfImage, InputQBallImageType, 0); + itkSetInputMacro(OdfImage, InputQBallImageType, 0) + itkGetInputMacro(OdfImage, InputQBallImageType, 0) /** Set/Get the Input mask image **/ - itkSetMacro(MaskImage, MaskImageTypePointer); - itkGetMacro(MaskImage, MaskImageTypePointer); + itkSetMacro(MaskImage, MaskImageTypePointer) + itkGetMacro(MaskImage, MaskImageTypePointer) - itkSetMacro(GfaImage, GfaImageTypePointer); - itkGetMacro(GfaImage, GfaImageTypePointer); + itkSetMacro(GfaImage, GfaImageTypePointer) + itkGetMacro(GfaImage, GfaImageTypePointer) - itkGetMacro(NumParticles, unsigned long); - itkGetMacro(NumConnections, unsigned long); - itkGetMacro(NumAcceptedFibers, int); - itkGetMacro(ProposalAcceptance, float); - itkGetMacro(Steps, unsigned int); + 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(); /** 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; }; virtual void Update(){ this->GenerateData(); } FiberPolyDataType GetFiberBundle(); float GetMemoryUsage(); bool EstimateParticleWeight(); protected: GibbsTrackingFilter(); virtual ~GibbsTrackingFilter(); void ComputeFiberCorrelation(); void ComputeFiberCorrelationOriginal(); void BuildFibers(float* points, int numPoints); // Input Images typename InputQBallImageType::Pointer m_ItkQBallImage; typename MaskImageType::Pointer m_MaskImage; typename GfaImageType::Pointer m_GfaImage; // 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 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; - RJMCMC* m_Sampler; + 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/Modules/DiffusionImaging/files.cmake b/Modules/DiffusionImaging/files.cmake index 94f64e08b7..df82fc2d51 100644 --- a/Modules/DiffusionImaging/files.cmake +++ b/Modules/DiffusionImaging/files.cmake @@ -1,229 +1,238 @@ 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 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 )