diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkFiberBuilder.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkFiberBuilder.cpp index f5add08263..e9d3fdfb7b 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkFiberBuilder.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkFiberBuilder.cpp @@ -1,138 +1,138 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkFiberBuilder.h" using namespace mitk; FiberBuilder::FiberBuilder(ParticleGrid* grid, ItkFloatImageType* image) { m_Grid = grid; m_Image = image; m_FiberLength = 0; } FiberBuilder::~FiberBuilder() { } vtkSmartPointer FiberBuilder::iterate(int minFiberLength) { m_VtkCellArray = vtkSmartPointer::New(); m_VtkPoints = vtkSmartPointer::New(); int cur_label = 1; int numFibers = 0; m_FiberLength = 0; for (int k = 0; k < m_Grid->m_NumParticles;k++) { Particle *dp = m_Grid->GetParticle(k); if (dp->label == 0) { vtkSmartPointer container = vtkSmartPointer::New(); dp->label = cur_label; LabelPredecessors(dp, -1, container); LabelSuccessors(dp, 1, container); cur_label++; if(m_FiberLength >= minFiberLength) { m_VtkCellArray->InsertNextCell(container); numFibers++; } m_FiberLength = 0; } } for (int k = 0; k < m_Grid->m_NumParticles;k++) { Particle *dp = m_Grid->GetParticle(k); dp->label = 0; } vtkSmartPointer fiberPolyData = vtkSmartPointer::New(); fiberPolyData->SetPoints(m_VtkPoints); fiberPolyData->SetLines(m_VtkCellArray); return fiberPolyData; } void FiberBuilder::LabelPredecessors(Particle* p, int ep, vtkPolyLine* container) { Particle* p2 = NULL; if (ep==1) p2 = m_Grid->GetParticle(p->pID); else p2 = m_Grid->GetParticle(p->mID); if (p2!=NULL && p2->label==0) { p2->label = 1; // assign particle to current fiber if (p2->pID==p->ID) LabelPredecessors(p2, -1, container); else if (p2->mID==p->ID) LabelPredecessors(p2, 1, container); else std::cout << "FiberBuilder: connection inconsistent (LabelPredecessors)" << std::endl; } AddPoint(p, container); } void FiberBuilder::LabelSuccessors(Particle* p, int ep, vtkPolyLine* container) { AddPoint(p, container); Particle* p2 = NULL; if (ep==1) p2 = m_Grid->GetParticle(p->pID); else p2 = m_Grid->GetParticle(p->mID); if (p2!=NULL && p2->label==0) { p2->label = 1; // assign particle to current fiber if (p2->pID==p->ID) LabelSuccessors(p2, -1, container); else if (p2->mID==p->ID) LabelSuccessors(p2, 1, container); else std::cout << "FiberBuilder: connection inconsistent (LabelPredecessors)" << std::endl; } } void FiberBuilder::AddPoint(Particle *dp, vtkSmartPointer container) { if (dp->label!=1) return; dp->label = 2; itk::ContinuousIndex index; - index[0] = dp->pos[0]/m_Image->GetSpacing()[0]-0.5; - index[1] = dp->pos[1]/m_Image->GetSpacing()[1]-0.5; - index[2] = dp->pos[2]/m_Image->GetSpacing()[2]-0.5; + index[0] = dp->GetPos()[0]/m_Image->GetSpacing()[0]-0.5; + index[1] = dp->GetPos()[1]/m_Image->GetSpacing()[1]-0.5; + index[2] = dp->GetPos()[2]/m_Image->GetSpacing()[2]-0.5; itk::Point point; m_Image->TransformContinuousIndexToPhysicalPoint( index, point ); vtkIdType id = m_VtkPoints->InsertNextPoint(point.GetDataPointer()); container->GetPointIds()->InsertNextId(id); if(container->GetNumberOfPoints()>1) m_FiberLength += m_LastPoint.EuclideanDistanceTo(point); m_LastPoint = point; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp index 0871309f2d..00790910ba 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkGibbsEnergyComputer.cpp @@ -1,218 +1,218 @@ #include #include #include #include #include #include using namespace mitk; GibbsEnergyComputer::GibbsEnergyComputer(ItkQBallImgType* qballImage, ItkFloatImageType* mask, ParticleGrid* particleGrid, SphereInterpolator* interpolator, ItkRandGenType* randGen) :EnergyComputer(mask, particleGrid, interpolator, randGen) { m_Image = qballImage; } GibbsEnergyComputer::~GibbsEnergyComputer() { } float GibbsEnergyComputer::EvaluateOdf(vnl_vector_fixed& pos, vnl_vector_fixed dir) { const int sampleSteps = 10; // evaluate ODF at 2*sampleSteps+1 positions along dir vnl_vector_fixed samplePos; // current position to evaluate float result = 0; // average of sampled ODF values int xint, yint, zint; // voxel containing samplePos // rotate particle direction according to image rotation dir = m_RotationMatrix*dir; // get interpolation for rotated direction m_SphereInterpolator->getInterpolation(dir); // sample ODF values along particle direction for (int i=-sampleSteps; i <= sampleSteps;i++) { samplePos = pos + (dir * m_ParticleLength) * ((float)i/sampleSteps); if (!m_UseTrilinearInterpolation) // image has not enough slices to use trilinear interpolation { ItkQBallImgType::IndexType index; index[0] = floor(pos[0]/m_Spacing[0]); index[1] = floor(pos[1]/m_Spacing[1]); index[2] = floor(pos[2]/m_Spacing[2]); if (m_Image->GetLargestPossibleRegion().IsInside(index)) { result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2]); } } else // use trilinear interpolation { float Rx = samplePos[0]/m_Spacing[0]-0.5; float Ry = samplePos[1]/m_Spacing[1]-0.5; float Rz = samplePos[2]/m_Spacing[2]-0.5; xint = floor(Rx); yint = floor(Ry); zint = floor(Rz); if (xint >= 0 && xint < m_Size[0]-1 && yint >= 0 && yint < m_Size[1]-1 && zint >= 0 && zint < m_Size[2]-1) { float xfrac = Rx-xint; float yfrac = Ry-yint; float zfrac = Rz-zint; ItkQBallImgType::IndexType index; float weight; weight = (1-xfrac)*(1-yfrac)*(1-zfrac); index[0] = xint; index[1] = yint; index[2] = zint; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(1-zfrac); index[0] = xint+1; index[1] = yint; index[2] = zint; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(1-zfrac); index[0] = xint; index[1] = yint+1; index[2] = zint; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(1-yfrac)*(zfrac); index[0] = xint; index[1] = yint; index[2] = zint+1; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(yfrac)*(1-zfrac); index[0] = xint+1; index[1] = yint+1; index[2] = zint; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (1-xfrac)*(yfrac)*(zfrac); index[0] = xint; index[1] = yint+1; index[2] = zint+1; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(1-yfrac)*(zfrac); index[0] = xint+1; index[1] = yint; index[2] = zint+1; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; weight = (xfrac)*(yfrac)*(zfrac); index[0] = xint+1; index[1] = yint+1; index[2] = zint+1; result += (m_Image->GetPixel(index)[m_SphereInterpolator->idx[0]-1]*m_SphereInterpolator->interpw[0] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[1]-1]*m_SphereInterpolator->interpw[1] + m_Image->GetPixel(index)[m_SphereInterpolator->idx[2]-1]* m_SphereInterpolator->interpw[2])*weight; } } } result /= (2*sampleSteps+1); // average result over taken samples return result; } float GibbsEnergyComputer::ComputeExternalEnergy(vnl_vector_fixed &R, vnl_vector_fixed &N, Particle *dp) { if (SpatProb(R) == 0) // check if position is inside mask return itk::NumericTraits::NonpositiveMin(); float odfVal = EvaluateOdf(R, N); // evaluate ODF in given direction float modelVal = 0; m_ParticleGrid->ComputeNeighbors(R); // retrieve neighbouring particles from particle grid Particle* neighbour = m_ParticleGrid->GetNextNeighbor(); while (neighbour!=NULL) // iterate over nieghbouring particles { if (dp != neighbour) // don't evaluate against itself { // see Reisert et al. "Global Reconstruction of Neuronal Fibers", MICCAI 2009 - float dot = fabs(dot_product(N,neighbour->dir)); + float dot = fabs(dot_product(N,neighbour->GetDir())); float bw = mbesseli0(dot); - float dpos = (neighbour->pos-R).squared_magnitude(); + float dpos = (neighbour->GetPos()-R).squared_magnitude(); float w = mexp(dpos*gamma_s); modelVal += w*(bw+m_ParticleChemicalPotential); w = mexp(dpos*gamma_reg_s); } neighbour = m_ParticleGrid->GetNextNeighbor(); } float energy = 2*(odfVal/m_ParticleWeight-modelVal) - (mbesseli0(1.0)+m_ParticleChemicalPotential); return energy*m_ExtStrength; } float GibbsEnergyComputer::ComputeInternalEnergy(Particle *dp) { float energy = 0; if (dp->pID != -1) // has predecessor energy += ComputeInternalEnergyConnection(dp,+1); if (dp->mID != -1) // has successor energy += ComputeInternalEnergyConnection(dp,-1); return energy; } float GibbsEnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1) { Particle *p2 = 0; int ep2 = 0; if (ep1 == 1) p2 = m_ParticleGrid->GetParticle(p1->pID); // get predecessor else p2 = m_ParticleGrid->GetParticle(p1->mID); // get successor // check in which direction the connected particle is pointing if (p2->mID == p1->ID) ep2 = -1; else if (p2->pID == p1->ID) ep2 = 1; else std::cout << "EnergyComputer: Connections are inconsistent!" << std::endl; return ComputeInternalEnergyConnection(p1,ep1,p2,ep2); } float GibbsEnergyComputer::ComputeInternalEnergyConnection(Particle *p1,int ep1, Particle *p2, int ep2) { // see Reisert et al. "Global Reconstruction of Neuronal Fibers", MICCAI 2009 - if ((dot_product(p1->dir,p2->dir))*ep1*ep2 > -m_CurvatureThreshold) // angle between particles is too sharp + if ((dot_product(p1->GetDir(),p2->GetDir()))*ep1*ep2 > -m_CurvatureThreshold) // angle between particles is too sharp return itk::NumericTraits::NonpositiveMin(); // calculate the endpoints of the two particles - vnl_vector_fixed endPoint1 = p1->pos + (p1->dir * (m_ParticleLength * ep1)); - vnl_vector_fixed endPoint2 = p2->pos + (p2->dir * (m_ParticleLength * ep2)); + vnl_vector_fixed endPoint1 = p1->GetPos() + (p1->GetDir() * (m_ParticleLength * ep1)); + vnl_vector_fixed endPoint2 = p2->GetPos() + (p2->GetDir() * (m_ParticleLength * ep2)); // check if endpoints are too far apart to connect if ((endPoint1-endPoint2).squared_magnitude() > m_SquaredParticleLength) return itk::NumericTraits::NonpositiveMin(); // calculate center point of the two particles - vnl_vector_fixed R = (p2->pos + p1->pos); R *= 0.5; + vnl_vector_fixed R = (p2->GetPos() + p1->GetPos()); R *= 0.5; // they are not allowed to connect if the mask image does not allow it if (SpatProb(R) == 0) return itk::NumericTraits::NonpositiveMin(); // get distances of endpoints to center point float norm1 = (endPoint1-R).squared_magnitude(); float norm2 = (endPoint2-R).squared_magnitude(); // calculate actual internal energy float energy = (m_ConnectionPotential-norm1-norm2)*m_IntStrength; return energy; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp index 4f0bf556f4..40b2058dfe 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkMetropolisHastingsSampler.cpp @@ -1,478 +1,478 @@ /*=================================================================== 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, EnergyComputer* enComp, ItkRandGenType* randGen, float curvThres) : m_ExTemp(0.01) , m_BirthProb(0.25) , m_DeathProb(0.05) , m_ShiftProb(0.15) , m_OptShiftProb(0.1) , m_ConnectionProb(0.45) , m_TractProb(0.5) , m_DelProb(0.1) , m_ChempotParticle(0.0) , m_AcceptedProposals(0) { m_RandGen = randGen; m_ParticleGrid = grid; m_EnergyComputer = enComp; m_ParticleLength = m_ParticleGrid->m_ParticleLength; m_DistanceThreshold = m_ParticleLength*m_ParticleLength; m_Sigma = m_ParticleLength/8.0; m_Gamma = 1/(m_Sigma*m_Sigma*2); m_Z = pow(2*M_PI*m_Sigma,3.0/2.0)*(M_PI*m_Sigma/m_ParticleLength); m_CurvatureThreshold = curvThres; m_StopProb = exp(-1/m_TractProb); } void MetropolisHastingsSampler::SetProbabilities(float birth, float death, float shift, float optShift, float connect) { m_BirthProb = birth; m_DeathProb = death; m_ShiftProb = shift; m_OptShiftProb = optShift; m_ConnectionProb = connect; float sum = m_BirthProb+m_DeathProb+m_ShiftProb+m_OptShiftProb+m_ConnectionProb; if (sum!=1 && sum>mitk::eps) { m_BirthProb /= sum; m_DeathProb /= sum; m_ShiftProb /= sum; m_OptShiftProb /= sum; m_ConnectionProb /= sum; } std::cout << "Update proposal probabilities:" << std::endl; std::cout << "Birth: " << m_BirthProb << std::endl; std::cout << "Death: " << m_DeathProb << std::endl; std::cout << "Shift: " << m_ShiftProb << std::endl; std::cout << "Optimal shift: " << m_OptShiftProb << std::endl; std::cout << "Connection: " << m_ConnectionProb << std::endl; } // update temperature of simulated annealing process void MetropolisHastingsSampler::SetTemperature(float val) { m_InTemp = val; m_Density = exp(-m_ChempotParticle/m_InTemp); } // add small random number drawn from gaussian to each vector element void MetropolisHastingsSampler::DistortVector(float sigma, vnl_vector_fixed& vec) { vec[0] += m_RandGen->GetNormalVariate(0.0, sigma); vec[1] += m_RandGen->GetNormalVariate(0.0, sigma); vec[2] += m_RandGen->GetNormalVariate(0.0, sigma); } // generate normalized random vector vnl_vector_fixed MetropolisHastingsSampler::GetRandomDirection() { vnl_vector_fixed vec; vec[0] = m_RandGen->GetNormalVariate(); vec[1] = m_RandGen->GetNormalVariate(); vec[2] = m_RandGen->GetNormalVariate(); vec.normalize(); return vec; } // generate actual proposal (birth, death, shift and connection of particle) void MetropolisHastingsSampler::MakeProposal() { float randnum = m_RandGen->GetVariate(); // Birth Proposal if (randnum < m_BirthProb) { vnl_vector_fixed R; m_EnergyComputer->DrawRandomPosition(R); vnl_vector_fixed N = GetRandomDirection(); Particle prop; - prop.pos = R; - prop.dir = N; + prop.GetPos() = R; + prop.GetDir() = N; float prob = m_Density * m_DeathProb /((m_BirthProb)*(m_ParticleGrid->m_NumParticles+1)); float ex_energy = m_EnergyComputer->ComputeExternalEnergy(R,N,0); float in_energy = m_EnergyComputer->ComputeInternalEnergy(&prop); prob *= exp((in_energy/m_InTemp+ex_energy/m_ExTemp)) ; if (prob > 1 || m_RandGen->GetVariate() < prob) { Particle *p = m_ParticleGrid->NewParticle(R); if (p!=0) { - p->pos = R; - p->dir = N; + p->GetPos() = R; + p->GetDir() = N; m_AcceptedProposals++; } } } // Death Proposal else if (randnum < m_BirthProb+m_DeathProb) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = m_RandGen->GetIntegerVariate()%m_ParticleGrid->m_NumParticles; Particle *dp = m_ParticleGrid->GetParticle(pnum); if (dp->pID == -1 && dp->mID == -1) { - float ex_energy = m_EnergyComputer->ComputeExternalEnergy(dp->pos,dp->dir,dp); + float ex_energy = m_EnergyComputer->ComputeExternalEnergy(dp->GetPos(),dp->GetDir(),dp); float in_energy = m_EnergyComputer->ComputeInternalEnergy(dp); float prob = m_ParticleGrid->m_NumParticles * (m_BirthProb) /(m_Density*m_DeathProb); //*SpatProb(dp->R); prob *= exp(-(in_energy/m_InTemp+ex_energy/m_ExTemp)) ; if (prob > 1 || m_RandGen->GetVariate() < prob) { m_ParticleGrid->RemoveParticle(pnum); m_AcceptedProposals++; } } } } // Shift Proposal else if (randnum < m_BirthProb+m_DeathProb+m_ShiftProb) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = m_RandGen->GetIntegerVariate()%m_ParticleGrid->m_NumParticles; Particle *p = m_ParticleGrid->GetParticle(pnum); Particle prop_p = *p; - DistortVector(m_Sigma, prop_p.pos); - DistortVector(m_Sigma/(2*m_ParticleLength), prop_p.dir); - prop_p.dir.normalize(); + DistortVector(m_Sigma, prop_p.GetPos()); + DistortVector(m_Sigma/(2*m_ParticleLength), prop_p.GetDir()); + prop_p.GetDir().normalize(); - float ex_energy = m_EnergyComputer->ComputeExternalEnergy(prop_p.pos,prop_p.dir,p) - - m_EnergyComputer->ComputeExternalEnergy(p->pos,p->dir,p); + float ex_energy = m_EnergyComputer->ComputeExternalEnergy(prop_p.GetPos(),prop_p.GetDir(),p) + - m_EnergyComputer->ComputeExternalEnergy(p->GetPos(),p->GetDir(),p); float in_energy = m_EnergyComputer->ComputeInternalEnergy(&prop_p) - m_EnergyComputer->ComputeInternalEnergy(p); float prob = exp(ex_energy/m_ExTemp+in_energy/m_InTemp); if (m_RandGen->GetVariate() < prob) { - vnl_vector_fixed Rtmp = p->pos; - vnl_vector_fixed Ntmp = p->dir; - p->pos = prop_p.pos; - p->dir = prop_p.dir; + vnl_vector_fixed Rtmp = p->GetPos(); + vnl_vector_fixed Ntmp = p->GetDir(); + p->GetPos() = prop_p.GetPos(); + p->GetDir() = prop_p.GetDir(); if (!m_ParticleGrid->TryUpdateGrid(pnum)) { - p->pos = Rtmp; - p->dir = Ntmp; + p->GetPos() = Rtmp; + p->GetDir() = Ntmp; } m_AcceptedProposals++; } } } // Optimal Shift Proposal else if (randnum < m_BirthProb+m_DeathProb+m_ShiftProb+m_OptShiftProb) { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = m_RandGen->GetIntegerVariate()%m_ParticleGrid->m_NumParticles; Particle *p = m_ParticleGrid->GetParticle(pnum); bool no_proposal = false; Particle prop_p = *p; if (p->pID != -1 && p->mID != -1) { Particle *plus = m_ParticleGrid->GetParticle(p->pID); int ep_plus = (plus->pID == p->ID)? 1 : -1; Particle *minus = m_ParticleGrid->GetParticle(p->mID); int ep_minus = (minus->pID == p->ID)? 1 : -1; - prop_p.pos = (plus->pos + plus->dir * (m_ParticleLength * ep_plus) + minus->pos + minus->dir * (m_ParticleLength * ep_minus)); - prop_p.pos *= 0.5; - prop_p.dir = plus->pos - minus->pos; - prop_p.dir.normalize(); + prop_p.GetPos() = (plus->GetPos() + plus->GetDir() * (m_ParticleLength * ep_plus) + minus->GetPos() + minus->GetDir() * (m_ParticleLength * ep_minus)); + prop_p.GetPos() *= 0.5; + prop_p.GetDir() = plus->GetPos() - minus->GetPos(); + prop_p.GetDir().normalize(); } else if (p->pID != -1) { Particle *plus = m_ParticleGrid->GetParticle(p->pID); int ep_plus = (plus->pID == p->ID)? 1 : -1; - prop_p.pos = plus->pos + plus->dir * (m_ParticleLength * ep_plus * 2); - prop_p.dir = plus->dir; + prop_p.GetPos() = plus->GetPos() + plus->GetDir() * (m_ParticleLength * ep_plus * 2); + prop_p.GetDir() = plus->GetDir(); } else if (p->mID != -1) { Particle *minus = m_ParticleGrid->GetParticle(p->mID); int ep_minus = (minus->pID == p->ID)? 1 : -1; - prop_p.pos = minus->pos + minus->dir * (m_ParticleLength * ep_minus * 2); - prop_p.dir = minus->dir; + prop_p.GetPos() = minus->GetPos() + minus->GetDir() * (m_ParticleLength * ep_minus * 2); + prop_p.GetDir() = minus->GetDir(); } else no_proposal = true; if (!no_proposal) { - float cos = dot_product(prop_p.dir, p->dir); - float p_rev = exp(-((prop_p.pos-p->pos).squared_magnitude() + (1-cos*cos))*m_Gamma)/m_Z; + float cos = dot_product(prop_p.GetDir(), p->GetDir()); + float p_rev = exp(-((prop_p.GetPos()-p->GetPos()).squared_magnitude() + (1-cos*cos))*m_Gamma)/m_Z; - float ex_energy = m_EnergyComputer->ComputeExternalEnergy(prop_p.pos,prop_p.dir,p) - - m_EnergyComputer->ComputeExternalEnergy(p->pos,p->dir,p); + float ex_energy = m_EnergyComputer->ComputeExternalEnergy(prop_p.GetPos(),prop_p.GetDir(),p) + - m_EnergyComputer->ComputeExternalEnergy(p->GetPos(),p->GetDir(),p); float in_energy = m_EnergyComputer->ComputeInternalEnergy(&prop_p) - m_EnergyComputer->ComputeInternalEnergy(p); float prob = exp(ex_energy/m_ExTemp+in_energy/m_InTemp)*m_ShiftProb*p_rev/(m_OptShiftProb+m_ShiftProb*p_rev); if (m_RandGen->GetVariate() < prob) { - vnl_vector_fixed Rtmp = p->pos; - vnl_vector_fixed Ntmp = p->dir; - p->pos = prop_p.pos; - p->dir = prop_p.dir; + vnl_vector_fixed Rtmp = p->GetPos(); + vnl_vector_fixed Ntmp = p->GetDir(); + p->GetPos() = prop_p.GetPos(); + p->GetDir() = prop_p.GetDir(); if (!m_ParticleGrid->TryUpdateGrid(pnum)) { - p->pos = Rtmp; - p->dir = Ntmp; + p->GetPos() = Rtmp; + p->GetDir() = Ntmp; } m_AcceptedProposals++; } } } } else { if (m_ParticleGrid->m_NumParticles > 0) { int pnum = m_RandGen->GetIntegerVariate()%m_ParticleGrid->m_NumParticles; Particle *p = m_ParticleGrid->GetParticle(pnum); EndPoint P; P.p = p; P.ep = (m_RandGen->GetVariate() > 0.5)? 1 : -1; RemoveAndSaveTrack(P); if (m_BackupTrack.m_Probability != 0) { MakeTrackProposal(P); float prob = (m_ProposalTrack.m_Energy-m_BackupTrack.m_Energy)/m_InTemp ; prob = exp(prob)*(m_BackupTrack.m_Probability * pow(m_DelProb,m_ProposalTrack.m_Length)) /(m_ProposalTrack.m_Probability * pow(m_DelProb,m_BackupTrack.m_Length)); if (m_RandGen->GetVariate() < prob) { ImplementTrack(m_ProposalTrack); m_AcceptedProposals++; } else { ImplementTrack(m_BackupTrack); } } else ImplementTrack(m_BackupTrack); } } } // establish connections between particles stored in input Track void MetropolisHastingsSampler::ImplementTrack(Track &T) { for (int k = 1; k < T.m_Length;k++) m_ParticleGrid->CreateConnection(T.track[k-1].p,T.track[k-1].ep,T.track[k].p,-T.track[k].ep); } // remove pending track from random particle, save it in m_BackupTrack and calculate its probability void MetropolisHastingsSampler::RemoveAndSaveTrack(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; m_BackupTrack.track[cnt] = Current; EndPoint Next; for (;;) { Next.p = 0; if (Current.ep == 1) { if (Current.p->pID != -1) { Next.p = m_ParticleGrid->GetParticle(Current.p->pID); Current.p->pID = -1; m_ParticleGrid->m_NumConnections--; } } else if (Current.ep == -1) { if (Current.p->mID != -1) { Next.p = m_ParticleGrid->GetParticle(Current.p->mID); Current.p->mID = -1; m_ParticleGrid->m_NumConnections--; } } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 3\n"); break; } if (Next.p == 0) // no successor { Next.ep = 0; // mark as empty successor break; } else { if (Next.p->pID == Current.p->ID) { Next.p->pID = -1; Next.ep = 1; } else if (Next.p->mID == Current.p->ID) { Next.p->mID = -1; Next.ep = -1; } else { fprintf(stderr,"MetropolisHastingsSampler_randshift: Connection inconsistent 4\n"); break; } } ComputeEndPointProposalDistribution(Current); AccumProb *= (m_SimpSamp.probFor(Next)); if (Next.p == 0) // no successor -> break break; energy += m_EnergyComputer->ComputeInternalEnergyConnection(Current.p,Current.ep,Next.p,Next.ep); Current = Next; Current.ep *= -1; cnt++; m_BackupTrack.track[cnt] = Current; if (m_RandGen->GetVariate() > m_DelProb) break; } m_BackupTrack.m_Energy = energy; m_BackupTrack.m_Probability = AccumProb; m_BackupTrack.m_Length = cnt+1; } // generate new track using kind of a local tracking starting from P in the given direction, store it in m_ProposalTrack and calculate its probability void MetropolisHastingsSampler::MakeTrackProposal(EndPoint P) { EndPoint Current = P; int cnt = 0; float energy = 0; float AccumProb = 1.0; m_ProposalTrack.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); int k = m_SimpSamp.draw(m_RandGen->GetVariate()); // stop tracking proposed if (k==0) break; EndPoint Next = m_SimpSamp.objs[k]; float probability = m_SimpSamp.probFor(k); // accumulate energy and proposal distribution energy += m_EnergyComputer->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 m_ProposalTrack.track[cnt++] = Current; } m_ProposalTrack.m_Energy = energy; m_ProposalTrack.m_Probability = AccumProb; m_ProposalTrack.m_Length = cnt; // clear labels for (int j = 0; j < m_ProposalTrack.m_Length;j++) m_ProposalTrack.track[j].p->label = 0; } // get neigbouring particles of P and calculate the according connection probabilities void MetropolisHastingsSampler::ComputeEndPointProposalDistribution(EndPoint P) { Particle *p = P.p; int ep = P.ep; float dist,dot; - vnl_vector_fixed R = p->pos + (p->dir * (ep*m_ParticleLength) ); + vnl_vector_fixed R = p->GetPos() + (p->GetDir() * (ep*m_ParticleLength) ); m_ParticleGrid->ComputeNeighbors(R); m_SimpSamp.clear(); m_SimpSamp.add(m_StopProb,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->pos - p2->dir * m_ParticleLength - R).squared_magnitude(); + dist = (p2->GetPos() - p2->GetDir() * m_ParticleLength - R).squared_magnitude(); if (dist < m_DistanceThreshold) { - dot = dot_product(p2->dir,p->dir) * ep; + dot = dot_product(p2->GetDir(),p->GetDir()) * ep; if (dot > m_CurvatureThreshold) { float en = m_EnergyComputer->ComputeInternalEnergyConnection(p,ep,p2,-1); m_SimpSamp.add(exp(en/m_TractProb),EndPoint(p2,-1)); } } } if (p2->pID == -1) { - dist = (p2->pos + p2->dir * m_ParticleLength - R).squared_magnitude(); + dist = (p2->GetPos() + p2->GetDir() * m_ParticleLength - R).squared_magnitude(); if (dist < m_DistanceThreshold) { - dot = dot_product(p2->dir,p->dir) * (-ep); + dot = dot_product(p2->GetDir(),p->GetDir()) * (-ep); if (dot > m_CurvatureThreshold) { float en = m_EnergyComputer->ComputeInternalEnergyConnection(p,ep,p2,+1); m_SimpSamp.add(exp(en/m_TractProb),EndPoint(p2,+1)); } } } } } } // return number of accepted proposals int MetropolisHastingsSampler::GetNumAcceptedProposals() { return m_AcceptedProposals; } diff --git a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h index a2e41268ef..03cb4c09db 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticle.h @@ -1,80 +1,94 @@ /*=================================================================== 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 { /** * \brief A particle is the basic element of the Gibbs fiber tractography method. */ class FiberTracking_EXPORT Particle { public: Particle() { label = 0; pID = -1; mID = -1; } ~Particle() { } + int gridindex; // index in the grid where it is living + int ID; // particle ID + int pID; // successor ID + int mID; // predecessor ID + unsigned char label; // label used in the fiber building process + + vnl_vector_fixed& GetPos() + { + return pos; + } + + vnl_vector_fixed& GetDir() + { + return dir; + } + +private: #pragma warning(push) #pragma warning(disable: 4251) + // this pragma ignores the following warning: + // warning C4251: 'mitk::Particle::pos' : class 'ATL::CStringT' needs to have dll-interface to be used by clients of class 'Particle' vnl_vector_fixed pos; // particle position (world coordinates. corner based voxels. not accounted for image rotation. vnl_vector_fixed dir; // normalized direction vector #pragma warning(pop) - int gridindex; // index in the grid where it is living - int ID; // particle ID - int pID; // successor ID - int mID; // predecessor ID - unsigned char label; // label used in the fiber building process }; class FiberTracking_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/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp index 7cefecb955..e3ef96efe2 100644 --- a/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp +++ b/Modules/DiffusionImaging/FiberTracking/Algorithms/GibbsTracking/mitkParticleGrid.cpp @@ -1,420 +1,420 @@ /*=================================================================== 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 particleLength, int cellCapacity) { // initialize counters m_NumParticles = 0; m_NumConnections = 0; m_NumCellOverflows = 0; m_ParticleLength = particleLength; // define isotropic grid from voxel spacing and particle length float cellSize = 2*m_ParticleLength; m_GridSize[0] = image->GetLargestPossibleRegion().GetSize()[0]*image->GetSpacing()[0]/cellSize +1; m_GridSize[1] = image->GetLargestPossibleRegion().GetSize()[1]*image->GetSpacing()[1]/cellSize +1; m_GridSize[2] = image->GetLargestPossibleRegion().GetSize()[2]*image->GetSpacing()[2]/cellSize +1; m_GridScale[0] = 1/cellSize; m_GridScale[1] = 1/cellSize; m_GridScale[2] = 1/cellSize; m_CellCapacity = cellCapacity; // maximum number of particles per grid cell m_ContainerCapacity = 100000; // initial particle container capacity unsigned long numCells = m_GridSize[0]*m_GridSize[1]*m_GridSize[2]; // number of grid cells unsigned long gridSize = numCells*m_CellCapacity; if ( (unsigned long)itk::NumericTraits::max() R) { if (m_NumParticles >= m_ContainerCapacity) { if (!ReallocateGrid()) return NULL; } int xint = int(R[0]*m_GridScale[0]); if (xint < 0) return NULL; if (xint >= m_GridSize[0]) return NULL; int yint = int(R[1]*m_GridScale[1]); if (yint < 0) return NULL; if (yint >= m_GridSize[1]) return NULL; int zint = int(R[2]*m_GridScale[2]); if (zint < 0) return NULL; if (zint >= m_GridSize[2]) return NULL; int idx = xint + m_GridSize[0]*(yint + m_GridSize[1]*zint); if (m_OccupationCount[idx] < m_CellCapacity) { Particle *p = &(m_Particles[m_NumParticles]); - p->pos = R; + p->GetPos() = R; p->mID = -1; p->pID = -1; m_NumParticles++; p->gridindex = m_CellCapacity*idx + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return p; } else { m_NumCellOverflows++; return NULL; } } bool ParticleGrid::TryUpdateGrid(int k) { Particle* p = &(m_Particles[k]); - int xint = int(p->pos[0]*m_GridScale[0]); + int xint = int(p->GetPos()[0]*m_GridScale[0]); if (xint < 0) return false; if (xint >= m_GridSize[0]) return false; - int yint = int(p->pos[1]*m_GridScale[1]); + int yint = int(p->GetPos()[1]*m_GridScale[1]); if (yint < 0) return false; if (yint >= m_GridSize[1]) return false; - int zint = int(p->pos[2]*m_GridScale[2]); + int zint = int(p->GetPos()[2]*m_GridScale[2]); if (zint < 0) return false; if (zint >= m_GridSize[2]) return false; int idx = xint + m_GridSize[0]*(yint+ zint*m_GridSize[1]); int cellidx = p->gridindex/m_CellCapacity; if (idx != cellidx) // cell has changed { if (m_OccupationCount[idx] < m_CellCapacity) { // remove from old position in grid; int grdindex = p->gridindex; m_Grid[grdindex] = m_Grid[cellidx*m_CellCapacity + m_OccupationCount[cellidx]-1]; m_Grid[grdindex]->gridindex = grdindex; m_OccupationCount[cellidx]--; // insert at new position in grid p->gridindex = idx*m_CellCapacity + m_OccupationCount[idx]; m_Grid[p->gridindex] = p; m_OccupationCount[idx]++; return true; } else { m_NumCellOverflows++; return false; } } return true; } void ParticleGrid::RemoveParticle(int k) { Particle* p = &(m_Particles[k]); int gridIndex = p->gridindex; int cellIdx = gridIndex/m_CellCapacity; int idx = gridIndex%m_CellCapacity; // remove pending connections if (p->mID != -1) DestroyConnection(p,-1); if (p->pID != -1) DestroyConnection(p,+1); // remove from grid if (idx < m_OccupationCount[cellIdx]-1) { m_Grid[gridIndex] = m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1]; m_Grid[cellIdx*m_CellCapacity+m_OccupationCount[cellIdx]-1] = NULL; m_Grid[gridIndex]->gridindex = gridIndex; } m_OccupationCount[cellIdx]--; // remove from container if (k < m_NumParticles-1) { Particle* last = &m_Particles[m_NumParticles-1]; // last particle // update connections of last particle because its index is changing if (last->mID!=-1) { if ( m_Particles[last->mID].mID == m_NumParticles-1 ) m_Particles[last->mID].mID = k; else if ( m_Particles[last->mID].pID == m_NumParticles-1 ) m_Particles[last->mID].pID = k; } if (last->pID!=-1) { if ( m_Particles[last->pID].mID == m_NumParticles-1 ) m_Particles[last->pID].mID = k; else if ( m_Particles[last->pID].pID == m_NumParticles-1 ) m_Particles[last->pID].pID = k; } m_Particles[k] = m_Particles[m_NumParticles-1]; // move very last particle to empty slot m_Particles[m_NumParticles-1].ID = m_NumParticles-1; // update ID of removed particle to match the index m_Particles[k].ID = k; // update ID of moved particle m_Grid[m_Particles[k].gridindex] = &m_Particles[k]; // update address of moved particle } m_NumParticles--; } void ParticleGrid::ComputeNeighbors(vnl_vector_fixed &R) { float xfrac = R[0]*m_GridScale[0]; float yfrac = R[1]*m_GridScale[1]; float zfrac = R[2]*m_GridScale[2]; int xint = int(xfrac); int yint = int(yfrac); int zint = int(zfrac); int dx = -1; if (xfrac-xint > 0.5) dx = 1; if (xint <= 0) { xint = 0; dx = 1; } if (xint >= m_GridSize[0]-1) { xint = m_GridSize[0]-1; dx = -1; } if (m_GridSize[0] <= 1) { dx = 0; } // Necessary with 2d images (bug 15416) int dy = -1; if (yfrac-yint > 0.5) dy = 1; if (yint <= 0) {yint = 0; dy = 1; } if (yint >= m_GridSize[1]-1) {yint = m_GridSize[1]-1; dy = -1;} if (m_GridSize[1] <= 1) { dy = 0; } // Necessary with 2d images (bug 15416) int dz = -1; if (zfrac-zint > 0.5) dz = 1; if (zint <= 0) {zint = 0; dz = 1; } if (zint >= m_GridSize[2]-1) {zint = m_GridSize[2]-1; dz = -1;} if (m_GridSize[2] <= 1) { dz = 0; } // Necessary with 2d images (bug 15416) m_NeighbourTracker.cellidx[0] = xint + m_GridSize[0]*(yint+zint*m_GridSize[1]); m_NeighbourTracker.cellidx[1] = m_NeighbourTracker.cellidx[0] + dx; m_NeighbourTracker.cellidx[2] = m_NeighbourTracker.cellidx[1] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[3] = m_NeighbourTracker.cellidx[2] - dx; m_NeighbourTracker.cellidx[4] = m_NeighbourTracker.cellidx[0] + dz*m_GridSize[0]*m_GridSize[1]; m_NeighbourTracker.cellidx[5] = m_NeighbourTracker.cellidx[4] + dx; m_NeighbourTracker.cellidx[6] = m_NeighbourTracker.cellidx[5] + dy*m_GridSize[0]; m_NeighbourTracker.cellidx[7] = m_NeighbourTracker.cellidx[6] - dx; m_NeighbourTracker.cellidx_c[0] = m_CellCapacity*m_NeighbourTracker.cellidx[0]; m_NeighbourTracker.cellidx_c[1] = m_CellCapacity*m_NeighbourTracker.cellidx[1]; m_NeighbourTracker.cellidx_c[2] = m_CellCapacity*m_NeighbourTracker.cellidx[2]; m_NeighbourTracker.cellidx_c[3] = m_CellCapacity*m_NeighbourTracker.cellidx[3]; m_NeighbourTracker.cellidx_c[4] = m_CellCapacity*m_NeighbourTracker.cellidx[4]; m_NeighbourTracker.cellidx_c[5] = m_CellCapacity*m_NeighbourTracker.cellidx[5]; m_NeighbourTracker.cellidx_c[6] = m_CellCapacity*m_NeighbourTracker.cellidx[6]; m_NeighbourTracker.cellidx_c[7] = m_CellCapacity*m_NeighbourTracker.cellidx[7]; m_NeighbourTracker.cellcnt = 0; m_NeighbourTracker.pcnt = 0; } Particle* ParticleGrid::GetNextNeighbor() { if (m_NeighbourTracker.pcnt < m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]]) { return m_Grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt] + (m_NeighbourTracker.pcnt++)]; } else { for(;;) { m_NeighbourTracker.cellcnt++; if (m_NeighbourTracker.cellcnt >= 8) return 0; if (m_OccupationCount[m_NeighbourTracker.cellidx[m_NeighbourTracker.cellcnt]] > 0) break; } m_NeighbourTracker.pcnt = 1; return m_Grid[m_NeighbourTracker.cellidx_c[m_NeighbourTracker.cellcnt]]; } } void ParticleGrid::CreateConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = P2->ID; else P1->pID = P2->ID; if (ep2 == -1) P2->mID = P1->ID; else P2->pID = P1->ID; m_NumConnections++; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1, Particle *P2, int ep2) { if (ep1 == -1) P1->mID = -1; else P1->pID = -1; if (ep2 == -1) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } void ParticleGrid::DestroyConnection(Particle *P1,int ep1) { Particle *P2 = 0; if (ep1 == 1) { P2 = &m_Particles[P1->pID]; P1->pID = -1; } else { P2 = &m_Particles[P1->mID]; P1->mID = -1; } if (P2->mID == P1->ID) P2->mID = -1; else P2->pID = -1; m_NumConnections--; } bool ParticleGrid::CheckConsistency() { for (int i=0; iID != i) { std::cout << "Particle ID error!" << std::endl; return false; } if (p->mID!=-1) { Particle* p2 = &m_Particles[p->mID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } if (p->pID!=-1) { Particle* p2 = &m_Particles[p->pID]; if (p2->mID!=p->ID && p2->pID!=p->ID) { std::cout << "Connection inconsistent!" << std::endl; return false; } } } return true; }