diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp index 60c0bf17c3..cef6166e84 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkFiberBuilder.cpp @@ -1,142 +1,142 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBuilder.h" using namespace mitk; FiberBuilder::FiberBuilder(ParticleGrid* grid, ItkFloatImageType* image) { m_Grid = grid; particles = m_Grid->m_Particles; m_Image = image; m_FiberLength = 0; } FiberBuilder::~FiberBuilder() { } vtkSmartPointer FiberBuilder::iterate(int minFiberLength) { m_VtkCellArray = vtkSmartPointer::New(); m_VtkPoints = vtkSmartPointer::New(); int cur_label = 1; int numFibers = 0; m_FiberLength = 0; for (int k = 0; k < m_Grid->m_NumParticles;k++) { Particle *dp = &(particles[k]); if (dp->label == 0) { vtkSmartPointer container = vtkSmartPointer::New(); dp->label = cur_label; dp->numerator = 0; labelPredecessors(dp, container); labelSuccessors(dp, container); cur_label++; if(m_FiberLength >= minFiberLength) { m_VtkCellArray->InsertNextCell(container); numFibers++; } m_FiberLength = 0; } } for (int k = 0; k < m_Grid->m_NumParticles;k++) { Particle *dp = &(particles[k]); dp->inserted = false; dp->label = 0; } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(m_VtkPoints); fiberPolyData->SetLines(m_VtkCellArray); return fiberPolyData; } void FiberBuilder::AddPoint(Particle *dp, vtkSmartPointer container) { if (dp->inserted) return; dp->inserted = true; itk::ContinuousIndex index; index[0] = dp->R[0]/m_Image->GetSpacing()[0]; index[1] = dp->R[1]/m_Image->GetSpacing()[1]; index[2] = dp->R[2]/m_Image->GetSpacing()[2]; itk::Point point; m_Image->TransformContinuousIndexToPhysicalPoint( index, point ); vtkIdType id = m_VtkPoints->InsertNextPoint(point.GetDataPointer()); container->GetPointIds()->InsertNextId(id); if(container->GetNumberOfPoints()>1) m_FiberLength += m_LastPoint.EuclideanDistanceTo(point); m_LastPoint = point; } void FiberBuilder::labelPredecessors(Particle *dp, vtkSmartPointer container) { if (dp->mID != -1 && dp->mID!=dp->ID) { - if (dp->ID!=particles[m_Grid->Id2Index(dp->mID)].pID) + if (dp->ID!=m_Grid->GetParticle(dp->mID)->pID) { - if (dp->ID==particles[m_Grid->Id2Index(dp->mID)].mID) + if (dp->ID==m_Grid->GetParticle(dp->mID)->mID) { - int tmp = particles[m_Grid->Id2Index(dp->mID)].pID; - particles[m_Grid->Id2Index(dp->mID)].pID = particles[m_Grid->Id2Index(dp->mID)].mID; - particles[m_Grid->Id2Index(dp->mID)].mID = tmp; + int tmp = m_Grid->GetParticle(dp->mID)->pID; + m_Grid->GetParticle(dp->mID)->pID = m_Grid->GetParticle(dp->mID)->mID; + m_Grid->GetParticle(dp->mID)->mID = tmp; } } - if (particles[m_Grid->Id2Index(dp->mID)].label == 0) + if (m_Grid->GetParticle(dp->mID)->label == 0) { - particles[m_Grid->Id2Index(dp->mID)].label = dp->label; - particles[m_Grid->Id2Index(dp->mID)].numerator = dp->numerator-1; - labelPredecessors(&(particles[m_Grid->Id2Index(dp->mID)]), container); + m_Grid->GetParticle(dp->mID)->label = dp->label; + m_Grid->GetParticle(dp->mID)->numerator = dp->numerator-1; + labelPredecessors(m_Grid->GetParticle(dp->mID), container); } } AddPoint(dp, container); } void FiberBuilder::labelSuccessors(Particle *dp, vtkSmartPointer container) { AddPoint(dp, container); if (dp->pID != -1 && dp->pID!=dp->ID) { - if (dp->ID!=particles[m_Grid->Id2Index(dp->pID)].mID) + if (dp->ID!=m_Grid->GetParticle(dp->pID)->mID) { - if (dp->ID==particles[m_Grid->Id2Index(dp->pID)].pID) + if (dp->ID==m_Grid->GetParticle(dp->pID)->pID) { - int tmp = particles[m_Grid->Id2Index(dp->pID)].pID; - particles[m_Grid->Id2Index(dp->pID)].pID = particles[m_Grid->Id2Index(dp->pID)].mID; - particles[m_Grid->Id2Index(dp->pID)].mID = tmp; + int tmp = m_Grid->GetParticle(dp->pID)->pID; + m_Grid->GetParticle(dp->pID)->pID = m_Grid->GetParticle(dp->pID)->mID; + m_Grid->GetParticle(dp->pID)->mID = tmp; } } - if (particles[m_Grid->Id2Index(dp->pID)].label == 0) + if (m_Grid->GetParticle(dp->pID)->label == 0) { - particles[m_Grid->Id2Index(dp->pID)].label = dp->label; - particles[m_Grid->Id2Index(dp->pID)].numerator = dp->numerator+1; - labelSuccessors(&(particles[m_Grid->Id2Index(dp->pID)]), container); + m_Grid->GetParticle(dp->pID)->label = dp->label; + m_Grid->GetParticle(dp->pID)->numerator = dp->numerator+1; + labelSuccessors(m_Grid->GetParticle(dp->pID), container); } } } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp index 194a357213..99d9bd7c7b 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.cpp @@ -1,646 +1,628 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkMetropolisHastingsSampler.h" using namespace mitk; MetropolisHastingsSampler::MetropolisHastingsSampler(ParticleGrid* grid, MTRand* randGen) : m_NumAttributes(0) , m_AcceptedProposals(0) { mtrand = randGen; m_ParticleGrid = grid; m_Iterations = 0; m_AcceptedProposals = 0; externalEnergy = 0; internalEnergy = 0; } void MetropolisHastingsSampler::SetEnergyComputer(EnergyComputer *e) { enc = e; } void MetropolisHastingsSampler::Iterate(float* acceptance, unsigned long* numCon, unsigned long* numPart, bool *abort) { m_AcceptedProposals = 0; for (int it = 0; it < m_Iterations;it++) { if (*abort) break; IterateOneStep(); *numCon = m_ParticleGrid->m_NumConnections; *numPart = m_ParticleGrid->m_NumParticles; } *acceptance = (float)m_AcceptedProposals/m_Iterations; } -void MetropolisHastingsSampler::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->Id2Index(p->mID); - npoints[m_NumAttributes*k+9] = m_ParticleGrid->Id2Index(p->pID); - } -} - void MetropolisHastingsSampler::SetParameters(float Temp, int numit, float plen, float curv_hardthres, float chempot_particle) { m_Iterations = numit; p_birth = 0.25; p_death = 0.05; p_shift = 0.15; p_shiftopt = 0.1; p_con = 0.45; p_cap = 0.0; m_ChempotParticle = chempot_particle; float sum = p_birth+p_death+p_shift+p_shiftopt+p_con; p_birth /= sum; p_death /= sum; p_shift /= sum; p_shiftopt /= sum; T_in = Temp; T_ex = 0.01; dens = exp(-chempot_particle/T_in); len_def = plen; len_sig = 0.0; cap_def = 1.0; cap_sig = 0.0; // shift proposal sigma_g = len_def/8.0; gamma_g = 1/(sigma_g*sigma_g*2); Z_g = pow(2*M_PI*sigma_g,3.0/2.0)*(M_PI*sigma_g/len_def); // conn proposal dthres = len_def; nthres = curv_hardthres; T_prop = 0.5; dthres *= dthres; stopprobability = exp(-1/T_prop); del_prob = 0.1; } void MetropolisHastingsSampler::SetTemperature(float temp) { T_in = temp; dens = exp(-m_ChempotParticle/T_in); } vnl_vector_fixed MetropolisHastingsSampler::distortn(float sigma, vnl_vector_fixed& vec) { vec[0] += sigma*mtrand->frandn(); vec[1] += sigma*mtrand->frandn(); vec[2] += sigma*mtrand->frandn(); } vnl_vector_fixed MetropolisHastingsSampler::rand_sphere() { vnl_vector_fixed vec; vec[0] += mtrand->frandn(); vec[1] += mtrand->frandn(); vec[2] += mtrand->frandn(); vec.normalize(); return vec; } //vnl_vector_fixed MetropolisHastingsSampler::rand(float w, float h, float d) //{ // vnl_vector_fixed vec; // vec[0] = mtrand->frand()*w; // vec[1] = mtrand->frand()*h; // vec[2] = mtrand->frand()*d; // return vec; //} void MetropolisHastingsSampler::IterateOneStep() { float randnum = mtrand->frand(); //randnum = 0; /////////////////////////////////////////////////////////////// //////// Birth Proposal /////////////////////////////////////////////////////////////// if (randnum < p_birth) { #ifdef TIMING tic(&birthproposal_time); birthstats.propose(); #endif vnl_vector_fixed R; enc->drawSpatPosition(R); //fprintf(stderr,"drawn: %f, %f, %f\n",R[0],R[1],R[2]); //R.setXYZ(20.5*3, 35.5*3, 1.5*3); vnl_vector_fixed N = rand_sphere(); //N.setXYZ(1,0,0); float cap = cap_def - cap_sig*mtrand->frand(); float len = len_def;// + len_sig*(mtrand->frand()-0.5); Particle prop; prop.R = R; prop.N = N; prop.cap = cap; prop.len = len; float prob = dens * p_death /((p_birth)*(m_ParticleGrid->m_NumParticles+1)); float ex_energy = enc->computeExternalEnergy(R,N,cap,len,0); float in_energy = enc->computeInternalEnergy(&prop); prob *= exp((in_energy/T_in+ex_energy/T_ex)) ; if (prob > 1 || mtrand->frand() < prob) { Particle *p = m_ParticleGrid->newParticle(R); if (p!=0) { p->R = R; p->N = N; p->cap = cap; p->len = len; #ifdef TIMING birthstats.accepted(); #endif m_AcceptedProposals++; } } #ifdef TIMING toc(&birthproposal_time); #endif } /////////////////////////////////////////////////////////////// //////// Death Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death) { if (m_ParticleGrid->m_NumParticles > 0) { #ifdef TIMING tic(&deathproposal_time); deathstats.propose(); #endif int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *dp = &(m_ParticleGrid->m_Particles[pnum]); if (dp->pID == -1 && dp->mID == -1) { float ex_energy = enc->computeExternalEnergy(dp->R,dp->N,dp->cap,dp->len,dp); float in_energy = enc->computeInternalEnergy(dp); float prob = m_ParticleGrid->m_NumParticles * (p_birth) /(dens*p_death); //*SpatProb(dp->R); prob *= exp(-(in_energy/T_in+ex_energy/T_ex)) ; if (prob > 1 || mtrand->frand() < prob) { m_ParticleGrid->RemoveParticle(pnum); #ifdef TIMING deathstats.accepted(); #endif m_AcceptedProposals++; } } #ifdef TIMING toc(&deathproposal_time); #endif } } /////////////////////////////////////////////////////////////// //////// Cap change Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death+p_cap) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); Particle prop_p = *p; prop_p.cap = cap_def - cap_sig*mtrand->frand(); float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,prop_p.cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); //float in_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) // - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float prob = exp(ex_energy/T_ex); // * SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand->frand() < prob) { p->cap = prop_p.cap; m_AcceptedProposals++; } } } /////////////////////////////////////////////////////////////// //////// Shift Proposal /////////////////////////////////////////////////////////////// else if (randnum < p_birth+p_death+p_shift+p_cap) { float energy = 0; if (m_ParticleGrid->m_NumParticles > 0) { #ifdef TIMING tic(&shiftproposal_time); shiftstats.propose(); #endif int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); Particle prop_p = *p; distortn(sigma_g, prop_p.R); distortn(sigma_g/(2*p->len), prop_p.N); prop_p.N.normalize(); float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float in_energy = enc->computeInternalEnergy(&prop_p) - enc->computeInternalEnergy(p); float prob = exp(ex_energy/T_ex+in_energy/T_in); // * SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand->frand() < prob) { vnl_vector_fixed Rtmp = p->R; vnl_vector_fixed Ntmp = p->N; p->R = prop_p.R; p->N = prop_p.N; if (!m_ParticleGrid->TryUpdateGrid(pnum)) { p->R = Rtmp; p->N = Ntmp; } #ifdef TIMING shiftstats.accepted(); #endif m_AcceptedProposals++; } #ifdef TIMING toc(&shiftproposal_time); #endif } } else if (randnum < p_birth+p_death+p_shift+p_shiftopt+p_cap) { float energy = 0; if (m_ParticleGrid->m_NumParticles > 0) { int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); bool no_proposal = false; Particle prop_p = *p; if (p->pID != -1 && p->mID != -1) { Particle *plus = m_ParticleGrid->m_AddressContainer[p->pID]; int ep_plus = (plus->pID == p->ID)? 1 : -1; Particle *minus = m_ParticleGrid->m_AddressContainer[p->mID]; int ep_minus = (minus->pID == p->ID)? 1 : -1; prop_p.R = (plus->R + plus->N * (plus->len * ep_plus) + minus->R + minus->N * (minus->len * ep_minus)); prop_p.R *= 0.5; prop_p.N = plus->R - minus->R; prop_p.N.normalize(); } else if (p->pID != -1) { Particle *plus = m_ParticleGrid->m_AddressContainer[p->pID]; int ep_plus = (plus->pID == p->ID)? 1 : -1; prop_p.R = plus->R + plus->N * (plus->len * ep_plus * 2); prop_p.N = plus->N; } else if (p->mID != -1) { Particle *minus = m_ParticleGrid->m_AddressContainer[p->mID]; int ep_minus = (minus->pID == p->ID)? 1 : -1; prop_p.R = minus->R + minus->N * (minus->len * ep_minus * 2); prop_p.N = minus->N; } else no_proposal = true; if (!no_proposal) { float cos = dot_product(prop_p.N, p->N); float p_rev = exp(-((prop_p.R-p->R).squared_magnitude() + (1-cos*cos))*gamma_g)/Z_g; float ex_energy = enc->computeExternalEnergy(prop_p.R,prop_p.N,p->cap,p->len,p) - enc->computeExternalEnergy(p->R,p->N,p->cap,p->len,p); float in_energy = enc->computeInternalEnergy(&prop_p) - enc->computeInternalEnergy(p); float prob = exp(ex_energy/T_ex+in_energy/T_in)*p_shift*p_rev/(p_shiftopt+p_shift*p_rev); //* SpatProb(p->R) / SpatProb(prop_p.R); if (mtrand->frand() < prob) { vnl_vector_fixed Rtmp = p->R; vnl_vector_fixed Ntmp = p->N; p->R = prop_p.R; p->N = prop_p.N; if (!m_ParticleGrid->TryUpdateGrid(pnum)) { p->R = Rtmp; p->N = Ntmp; } m_AcceptedProposals++; } } } } else { if (m_ParticleGrid->m_NumParticles > 0) { #ifdef TIMING tic(&connproposal_time); connstats.propose(); #endif int pnum = rand()%m_ParticleGrid->m_NumParticles; Particle *p = &(m_ParticleGrid->m_Particles[pnum]); EndPoint P; P.p = p; P.ep = (mtrand->frand() > 0.5)? 1 : -1; RemoveAndSaveTrack(P); if (TrackBackup.proposal_probability != 0) { MakeTrackProposal(P); float prob = (TrackProposal.energy-TrackBackup.energy)/T_in ; // prob = exp(prob)*(TrackBackup.proposal_probability) // /(TrackProposal.proposal_probability); prob = exp(prob)*(TrackBackup.proposal_probability * pow(del_prob,TrackProposal.length)) /(TrackProposal.proposal_probability * pow(del_prob,TrackBackup.length)); if (mtrand->frand() < prob) { ImplementTrack(TrackProposal); #ifdef TIMING connstats.accepted(); #endif m_AcceptedProposals++; } else { ImplementTrack(TrackBackup); } } else ImplementTrack(TrackBackup); #ifdef TIMING toc(&connproposal_time); #endif } } } void MetropolisHastingsSampler::ImplementTrack(Track &T) { for (int k = 1; k < T.length;k++) { m_ParticleGrid->CreateConnection(T.track[k-1].p,T.track[k-1].ep,T.track[k].p,-T.track[k].ep); } } void MetropolisHastingsSampler::RemoveAndSaveTrack(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; TrackBackup.track[cnt] = Current; EndPoint Next; for (;;) { Next.p = 0; if (Current.ep == 1) { if (Current.p->pID != -1) { Next.p = m_ParticleGrid->m_AddressContainer[Current.p->pID]; Current.p->pID = -1; m_ParticleGrid->m_NumConnections--; } } else if (Current.ep == -1) { if (Current.p->mID != -1) { Next.p = m_ParticleGrid->m_AddressContainer[Current.p->mID]; Current.p->mID = -1; m_ParticleGrid->m_NumConnections--; } } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 3\n"); break; } if (Next.p == 0) // no successor { Next.ep = 0; // mark as empty successor break; } else { if (Next.p->pID == Current.p->ID) { Next.p->pID = -1; Next.ep = 1; } else if (Next.p->mID == Current.p->ID) { Next.p->mID = -1; Next.ep = -1; } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 4\n"); break; } } ComputeEndPointProposalDistribution(Current); AccumProb *= (simpsamp.probFor(Next)); if (Next.p == 0) // no successor -> break break; energy += enc->computeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); Current = Next; Current.ep *= -1; cnt++; TrackBackup.track[cnt] = Current; if (mtrand->rand() > del_prob) { break; } } TrackBackup.energy = energy; TrackBackup.proposal_probability = AccumProb; TrackBackup.length = cnt+1; } void MetropolisHastingsSampler::MakeTrackProposal(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; TrackProposal.track[cnt++] = Current; Current.p->label = 1; for (;;) { // next candidate is already connected if ((Current.ep == 1 && Current.p->pID != -1) || (Current.ep == -1 && Current.p->mID != -1)) break; // track too long if (cnt > 250) break; ComputeEndPointProposalDistribution(Current); // // no candidates anymore // if (simpsamp.isempty()) // break; int k = simpsamp.draw(mtrand->frand()); // stop tracking proposed if (k==0) break; EndPoint Next = simpsamp.objs[k]; float probability = simpsamp.probFor(k); // accumulate energy and proposal distribution energy += enc->computeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); AccumProb *= probability; // track to next endpoint Current = Next; Current.ep *= -1; Current.p->label = 1; // put label to avoid loops TrackProposal.track[cnt++] = Current; } TrackProposal.energy = energy; TrackProposal.proposal_probability = AccumProb; TrackProposal.length = cnt; // clear labels for (int j = 0; j < TrackProposal.length;j++) TrackProposal.track[j].p->label = 0; } void MetropolisHastingsSampler::ComputeEndPointProposalDistribution(EndPoint P) { Particle *p = P.p; int ep = P.ep; float dist,dot; vnl_vector_fixed R = p->R + (p->N * (ep*p->len) ); m_ParticleGrid->ComputeNeighbors(R); simpsamp.clear(); simpsamp.add(stopprobability,EndPoint(0,0)); for (;;) { Particle *p2 = m_ParticleGrid->GetNextNeighbor(); if (p2 == 0) break; if (p!=p2 && p2->label == 0) { if (p2->mID == -1) { dist = (p2->R - p2->N * p2->len - R).squared_magnitude(); if (dist < dthres) { dot = dot_product(p2->N,p->N) * ep; if (dot > nthres) { float en = enc->computeInternalEnergyConnection(p,ep,p2,-1); simpsamp.add(exp(en/T_prop),EndPoint(p2,-1)); } } } if (p2->pID == -1) { dist = (p2->R + p2->N * p2->len - R).squared_magnitude(); if (dist < dthres) { dot = dot_product(p2->N,p->N) * (-ep); if (dot > nthres) { float en = enc->computeInternalEnergyConnection(p,ep,p2,+1); simpsamp.add(exp(en/T_prop),EndPoint(p2,+1)); } } } } } } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h index af674b2e50..49de82fb54 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkMetropolisHastingsSampler.h @@ -1,113 +1,111 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _SAMPLER #define _SAMPLER // MITK #include #include "mitkParticleGrid.h" #include "mitkEnergyComputer.h" #include "MersenneTwister.h" #include "mitkSimpSamp.h" // ITK #include // MISC #include namespace mitk { class MitkDiffusionImaging_EXPORT MetropolisHastingsSampler { public: typedef itk::Image< float, 3 > ItkFloatImageType; ParticleGrid* m_ParticleGrid; const int* datasz; EnergyComputer* enc; int m_Iterations; float width; float height; float depth; int m_NumAttributes; int m_AcceptedProposals; float T_in ; float T_ex ; float dens; float p_birth; float p_death; float p_shift; float p_shiftopt; float p_cap; float p_con; 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(ParticleGrid* grid, MTRand* randGen); - void WriteOutParticles(float *npoints); - void SetEnergyComputer(EnergyComputer *e); void SetParameters(float Temp, int numit, float plen, float curv_hardthres, float chempot_particle); void SetTemperature(float temp); void Iterate(float* acceptance, unsigned long* numCon, unsigned long* numPart, bool *abort); void IterateOneStep(); void ImplementTrack(Track &T); void RemoveAndSaveTrack(EndPoint P); void MakeTrackProposal(EndPoint P); void ComputeEndPointProposalDistribution(EndPoint P); vnl_vector_fixed distortn(float sigma, vnl_vector_fixed& vec); vnl_vector_fixed rand_sphere(); // vnl_vector_fixed rand(float w, float h, float d); }; } #endif diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp index 1db8dcf8c0..04dcc9129f 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.cpp @@ -1,424 +1,429 @@ /*=================================================================== 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(ItkFloatImageType* image, float cellSize) { //// involving the container capacity = 0; m_Particles = 0; m_AddressContainer = 0; m_NumParticles = 0; m_NumConnections = 0; m_NumCellOverflows = 0; ////// involvin the grid nx = 0; ny = 0; nz = 0; csize = 0; gridsize = 0; grid = (Particle**) 0; occnt = (int*) 0; increase_step = 100000; m_Memory = 0; int width = image->GetLargestPossibleRegion().GetSize()[0]*image->GetSpacing()[0]; int height = image->GetLargestPossibleRegion().GetSize()[1]*image->GetSpacing()[1]; int depth = image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]; float cellcnt_x = (int)((float)width/cellSize) +1; float cellcnt_y = (int)((float)height/cellSize) +1; float cellcnt_z = (int)((float)depth/cellSize) +1; int cell_capacity = 512; //// involving the container capacity = 100000; m_Particles = (Particle*) malloc(sizeof(Particle)*capacity); m_AddressContainer = (Particle**) malloc(sizeof(Particle*)*capacity); if (m_Particles == 0 || m_AddressContainer == 0) { fprintf(stderr,"error: Out of Memory\n"); capacity = 0; } else { fprintf(stderr,"Allocated Memory for %i particles \n",capacity); } 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 } // involving the grid nx = cellcnt_x; ny = cellcnt_y; nz = cellcnt_z; csize = cell_capacity; gridsize = nx*ny*nz; m_Memory = (float)(sizeof(Particle*)*gridsize*csize)/1000000; grid = (Particle**) malloc(sizeof(Particle*)*gridsize*csize); occnt = (int*) malloc(sizeof(int)*gridsize); if (grid == 0 || occnt == 0) { fprintf(stderr,"error: Out of Memory\n"); capacity = 0; } for (i = 0;i < gridsize;i++) occnt[i] = 0; mulx = 1/cellSize; muly = 1/cellSize; mulz = 1/cellSize; } -int ParticleGrid::reallocate() +int ParticleGrid::ReallocateGrid() { 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::Id2Index(int ID) { if (ID == -1) return -1; else return (m_AddressContainer[ID] - m_Particles); } +Particle* ParticleGrid::GetParticle(int ID) +{ + return &m_Particles[Id2Index(ID)]; +} + 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) + if (ReallocateGrid() == -1) { fprintf(stderr,"out of Memory!!\n"); return 0; } } int xint = int(R[0]*mulx); if (xint < 0) { //fprintf(stderr,"error: out of grid\n"); return 0;} if (xint >= nx) { // fprintf(stderr,"error: out of grid\n"); return 0;} int yint = int(R[1]*muly); if (yint < 0) { //fprintf(stderr,"error: out of grid\n"); return 0;} if (yint >= ny) {// fprintf(stderr,"error: out of grid\n"); return 0;} int zint = int(R[2]*mulz); if (zint < 0) {// fprintf(stderr,"error: out of grid\n"); return 0;} if (zint >= nz) { //fprintf(stderr,"error: out of grid\n"); return 0;} int idx = xint + nx*(yint + ny*zint); if (occnt[idx] < csize) { Particle *p = &(m_Particles[m_NumParticles]); p->R = R; p->mID = -1; p->pID = -1; m_NumParticles++; p->gridindex = csize*idx + occnt[idx]; grid[p->gridindex] = p; occnt[idx]++; return p; } else { m_NumCellOverflows++; //fprintf(stderr,"error: cell overflow \n"); return 0; } } bool ParticleGrid::TryUpdateGrid(int k) { Particle* p = &(m_Particles[k]); /////// find new grid cell int xint = int(p->R[0]*mulx); if (xint < 0) { return false; } if (xint >= nx) { return false; } int yint = int(p->R[1]*muly); if (yint < 0) { return false; } if (yint >= ny) { return false; } int zint = int(p->R[2]*mulz); if (zint < 0) { return false; } if (zint >= nz) { return false; } int idx = xint + nx*(yint+ zint*ny); int cellidx = p->gridindex/csize; if (idx != cellidx) // cell has changed { if (occnt[idx] < csize) { // remove from old position in grid; int grdindex = p->gridindex; grid[grdindex] = grid[cellidx*csize + occnt[cellidx]-1]; grid[grdindex]->gridindex = grdindex; occnt[cellidx]--; // insert at new position in grid p->gridindex = idx*csize + occnt[idx]; grid[p->gridindex] = p; occnt[idx]++; return true; } else return false; } return true; } void ParticleGrid::RemoveParticle(int k) { Particle* p = &(m_Particles[k]); int grdindex = p->gridindex; int cellidx = grdindex/csize; int idx = grdindex%csize; // remove pending connections if (p->mID != -1) DestroyConnection(p,-1); if (p->pID != -1) DestroyConnection(p,+1); // remove from grid if (idx < occnt[cellidx]-1) { grid[grdindex] = grid[cellidx*csize + occnt[cellidx]-1]; grid[grdindex]->gridindex = grdindex; } occnt[cellidx]--; // remove from container if (kID; int move_ID = m_Particles[m_NumParticles-1].ID; *p = m_Particles[m_NumParticles-1]; // move very last particle to empty slot m_Particles[m_NumParticles-1].ID = todel_ID; // keep IDs unique grid[p->gridindex] = p; // keep gridindex consistent // permute address table m_AddressContainer[todel_ID] = &(m_Particles[m_NumParticles-1]); m_AddressContainer[move_ID] = p; } m_NumParticles--; } void ParticleGrid::ComputeNeighbors(vnl_vector_fixed &R) { float xfrac = R[0]*mulx; float yfrac = R[1]*muly; float zfrac = R[2]*mulz; int xint = int(xfrac); int yint = int(yfrac); int zint = int(zfrac); int dx = -1; if (xfrac-xint > 0.5) dx = 1; if (xint <= 0) { xint = 0; dx = 1; } if (xint >= nx-1) { xint = nx-1; dx = -1; } int dy = -1; if (yfrac-yint > 0.5) dy = 1; if (yint <= 0) {yint = 0; dy = 1; } if (yint >= ny-1) {yint = ny-1; dy = -1;} int dz = -1; if (zfrac-zint > 0.5) dz = 1; if (zint <= 0) {zint = 0; dz = 1; } if (zint >= nz-1) {zint = nz-1; dz = -1;} m_NeighbourTracker.cellidx[0] = xint + nx*(yint+zint*ny); m_NeighbourTracker.cellidx[1] = m_NeighbourTracker.cellidx[0] + dx; m_NeighbourTracker.cellidx[2] = m_NeighbourTracker.cellidx[1] + dy*nx; m_NeighbourTracker.cellidx[3] = m_NeighbourTracker.cellidx[2] - dx; m_NeighbourTracker.cellidx[4] = m_NeighbourTracker.cellidx[0] + dz*nx*ny; m_NeighbourTracker.cellidx[5] = m_NeighbourTracker.cellidx[4] + dx; m_NeighbourTracker.cellidx[6] = m_NeighbourTracker.cellidx[5] + dy*nx; m_NeighbourTracker.cellidx[7] = m_NeighbourTracker.cellidx[6] - dx; m_NeighbourTracker.cellidx_c[0] = csize*m_NeighbourTracker.cellidx[0]; m_NeighbourTracker.cellidx_c[1] = csize*m_NeighbourTracker.cellidx[1]; m_NeighbourTracker.cellidx_c[2] = csize*m_NeighbourTracker.cellidx[2]; m_NeighbourTracker.cellidx_c[3] = csize*m_NeighbourTracker.cellidx[3]; m_NeighbourTracker.cellidx_c[4] = csize*m_NeighbourTracker.cellidx[4]; m_NeighbourTracker.cellidx_c[5] = csize*m_NeighbourTracker.cellidx[5]; m_NeighbourTracker.cellidx_c[6] = csize*m_NeighbourTracker.cellidx[6]; m_NeighbourTracker.cellidx_c[7] = csize*m_NeighbourTracker.cellidx[7]; m_NeighbourTracker.cellcnt = 0; m_NeighbourTracker.pcnt = 0; } Particle* ParticleGrid::GetNextNeighbor() { if (m_NeighbourTracker.pcnt < occnt[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]]) { return grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt] + (m_NeighbourTracker.pcnt++)]; } else { for(;;) { m_NeighbourTracker.cellcnt++; if (m_NeighbourTracker.cellcnt >= 8) return 0; if (occnt[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]] > 0) break; } m_NeighbourTracker.pcnt = 1; return grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt]]; } } void ParticleGrid::CreateConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = P2->ID; else P1->pID = P2->ID; if (ep2 == -1) P2->mID = P1->ID; else P2->pID = P1->ID; m_NumConnections++; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = -1; else P1->pID = -1; if (ep2 == -1) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1) { Particle *P2 = 0; if (ep1 == 1) { P2 = m_AddressContainer[P1->pID]; P1->pID = -1; } else { P2 = m_AddressContainer[P1->mID]; P1->mID = -1; } if (P2->mID == P1->ID) { P2->mID = -1; } else { P2->pID = -1; } m_NumConnections--; } diff --git a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h index e69b74381d..d0c87b4f85 100644 --- a/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h +++ b/Modules/DiffusionImaging/Tractography/GibbsTracking/mitkParticleGrid.h @@ -1,123 +1,123 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _PARTICLEGRID #define _PARTICLEGRID // MITK #include "mitkParticle.h" #include "MitkDiffusionImagingExports.h" // ITK #include namespace mitk { class MitkDiffusionImaging_EXPORT ParticleGrid { - //////////////// Container public: typedef itk::Image< float, 3 > ItkFloatImageType; Particle* m_Particles; // particles in linear array Particle** m_AddressContainer; int m_NumParticles; // number of particles int m_NumConnections; // number of connections int m_NumCellOverflows; // number of cell overflows ParticleGrid(ItkFloatImageType* image, float cellSize); ~ParticleGrid(); - int reallocate(); - - int Id2Index(int ID); + Particle* GetParticle(int ID); Particle* newParticle(vnl_vector_fixed R); bool TryUpdateGrid(int k); void RemoveParticle(int k); void ComputeNeighbors(vnl_vector_fixed &R); Particle* GetNextNeighbor(); void CreateConnection(Particle *P1,int ep1, Particle *P2, int ep2); void DestroyConnection(Particle *P1,int ep1, Particle *P2, int ep2); void DestroyConnection(Particle *P1,int ep1); private: + int Id2Index(int ID); + int ReallocateGrid(); + 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