diff --git a/Modules/DiffusionImaging/Testing/mitkGibbsTrackingTest.cpp b/Modules/DiffusionImaging/Testing/mitkGibbsTrackingTest.cpp index 98a470b701..a17e1213d7 100644 --- a/Modules/DiffusionImaging/Testing/mitkGibbsTrackingTest.cpp +++ b/Modules/DiffusionImaging/Testing/mitkGibbsTrackingTest.cpp @@ -1,99 +1,99 @@ /*=================================================================== 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 #include #include #include #include #include #include using namespace mitk; /**Documentation * Test for gibbs tracking filter */ int mitkGibbsTrackingTest(int argc, char* argv[]) { MITK_TEST_BEGIN("mitkGibbsTrackingTest"); MITK_TEST_CONDITION_REQUIRED(argc>5,"check for input data") QBallImage::Pointer mitkQballImage; Image::Pointer mitkMaskImage; mitk::FiberBundleX::Pointer fib1; try{ MITK_INFO << "Q-Ball image: " << argv[1]; MITK_INFO << "Mask image: " << argv[2]; MITK_INFO << "Parameter file: " << argv[3]; MITK_INFO << "Lut path: " << argv[4]; MITK_INFO << "Reference bundle: " << argv[5]; RegisterDiffusionImagingObjectFactory(); // test if fib1 can be read const std::string s1="", s2=""; std::vector infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[1], s1, s2, false ); mitkQballImage = dynamic_cast(infile.at(0).GetPointer()); MITK_TEST_CONDITION_REQUIRED(mitkQballImage.IsNotNull(),"check qball image") infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[2], s1, s2, false ); mitkMaskImage = dynamic_cast(infile.at(0).GetPointer()); MITK_TEST_CONDITION_REQUIRED(mitkMaskImage.IsNotNull(),"check mask image") infile = mitk::BaseDataIO::LoadBaseDataFromFile( argv[5], s1, s2, false ); fib1 = dynamic_cast(infile.at(0).GetPointer()); MITK_TEST_CONDITION_REQUIRED(fib1.IsNotNull(),"check fiber bundle") typedef itk::Vector OdfVectorType; typedef itk::Image OdfVectorImgType; typedef itk::Image MaskImgType; typedef itk::GibbsTrackingFilter GibbsTrackingFilterType; OdfVectorImgType::Pointer itk_qbi = OdfVectorImgType::New(); mitk::CastToItkImage(mitkQballImage, itk_qbi); MaskImgType::Pointer itk_mask = MaskImgType::New(); mitk::CastToItkImage(mitkMaskImage, itk_mask); GibbsTrackingFilterType::Pointer gibbsTracker = GibbsTrackingFilterType::New(); gibbsTracker->SetQBallImage(itk_qbi.GetPointer()); gibbsTracker->SetMaskImage(itk_mask); gibbsTracker->SetDuplicateImage(false); gibbsTracker->SetRandomSeed(1); - gibbsTracker->SetParameterFile(argv[3]); + gibbsTracker->SetLoadParameterFile(argv[3]); gibbsTracker->SetLutPath(argv[4]); gibbsTracker->Update(); mitk::FiberBundleX::Pointer fib2 = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); MITK_TEST_CONDITION_REQUIRED(fib1->Equals(fib2), "check if gibbs tracking has changed"); gibbsTracker->SetRandomSeed(0); gibbsTracker->Update(); fib2 = mitk::FiberBundleX::New(gibbsTracker->GetFiberBundle()); MITK_TEST_CONDITION_REQUIRED(!fib1->Equals(fib2), "check if gibbs tracking has changed after wrong seed"); } catch(...) { return EXIT_FAILURE; } // always end with this! MITK_TEST_END(); } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp index 95c38da0e0..5597642336 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.cpp @@ -1,416 +1,477 @@ /*=================================================================== 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" // MITK #include #include #include #include #include //#include #include #include // ITK #include #include +#include // MISC #include #include #include #include namespace itk{ template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::GibbsTrackingFilter(): m_StartTemperature(0.1), m_EndTemperature(0.001), m_Iterations(500000), m_ParticleWeight(0), m_ParticleWidth(0), m_ParticleLength(0), m_ConnectionPotential(10), m_InexBalance(0), m_ParticlePotential(0.2), m_MinFiberLength(10), m_AbortTracking(false), m_NumConnections(0), m_NumParticles(0), m_NumAcceptedFibers(0), m_CurrentStep(0), m_BuildFibers(false), m_Steps(10), m_ProposalAcceptance(0), m_CurvatureThreshold(0.7), m_DuplicateImage(true), m_RandomSeed(-1), - m_ParameterFile(""), + m_LoadParameterFile(""), m_LutPath("") { } template< class ItkQBallImageType > GibbsTrackingFilter< ItkQBallImageType >::~GibbsTrackingFilter() { } // fill output fiber bundle datastructure template< class ItkQBallImageType > typename GibbsTrackingFilter< ItkQBallImageType >::FiberPolyDataType GibbsTrackingFilter< ItkQBallImageType >::GetFiberBundle() { if (!m_AbortTracking) { m_BuildFibers = true; while (m_BuildFibers){} } return m_FiberPolyData; } template< class ItkQBallImageType > void GibbsTrackingFilter< ItkQBallImageType > ::EstimateParticleWeight() { MITK_INFO << "GibbsTrackingFilter: estimating particle weight"; float minSpacing; if(m_QBallImage->GetSpacing()[0]GetSpacing()[1] && m_QBallImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[0]; else if (m_QBallImage->GetSpacing()[1] < m_QBallImage->GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[1]; else minSpacing = m_QBallImage->GetSpacing()[2]; float m_ParticleLength = 1.5*minSpacing; float m_ParticleWidth = 0.5*minSpacing; // seed random generators Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = Statistics::MersenneTwisterRandomVariateGenerator::New(); if (m_RandomSeed>-1) randGen->SetSeed(m_RandomSeed); else randGen->SetSeed(); // instantiate all necessary components SphereInterpolator* interpolator = new SphereInterpolator(m_LutPath); ParticleGrid* particleGrid = new ParticleGrid(m_MaskImage, m_ParticleLength); GibbsEnergyComputer* encomp = new GibbsEnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); // EnergyComputer* encomp = new EnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); MetropolisHastingsSampler* sampler = new MetropolisHastingsSampler(particleGrid, encomp, randGen, m_CurvatureThreshold); float alpha = log(m_EndTemperature/m_StartTemperature); m_ParticleWeight = 0.01; int ppv = 0; // main loop int neededParts = 3000; while (ppvSetParameters(m_ParticleWeight,m_ParticleWidth,m_ConnectionPotential*m_ParticleLength*m_ParticleLength,m_CurvatureThreshold,m_InexBalance,m_ParticlePotential); for( int step = 0; step < 10; step++ ) { // update temperatur for simulated annealing process float temperature = m_StartTemperature * exp(alpha*(((1.0)*step)/((1.0)*10))); sampler->SetTemperature(temperature); for (unsigned long i=0; i<10000; i++) sampler->MakeProposal(); } ppv = particleGrid->m_NumParticles; particleGrid->ResetGrid(); } delete sampler; delete encomp; delete particleGrid; delete interpolator; MITK_INFO << "GibbsTrackingFilter: finished estimating particle weight"; } // perform global tracking template< class ItkQBallImageType > void GibbsTrackingFilter< ItkQBallImageType >::GenerateData() { + TimeProbe preClock; preClock.Start(); // check if input is qball or tensor image and generate qball if necessary if (m_QBallImage.IsNull() && m_TensorImage.IsNotNull()) { TensorImageToQBallImageFilter::Pointer filter = TensorImageToQBallImageFilter::New(); filter->SetInput( m_TensorImage ); filter->Update(); m_QBallImage = filter->GetOutput(); } else if (m_DuplicateImage) // generate local working copy of QBall image (if not disabled) { typedef itk::ImageDuplicator< ItkQBallImageType > DuplicateFilterType; typename DuplicateFilterType::Pointer duplicator = DuplicateFilterType::New(); duplicator->SetInputImage( m_QBallImage ); duplicator->Update(); m_QBallImage = duplicator->GetOutput(); } // perform mean subtraction on odfs typedef ImageRegionIterator< ItkQBallImageType > InputIteratorType; InputIteratorType it(m_QBallImage, m_QBallImage->GetLargestPossibleRegion() ); it.GoToBegin(); while (!it.IsAtEnd()) { itk::OrientationDistributionFunction odf(it.Get().GetDataPointer()); float mean = odf.GetMeanValue(); odf -= mean; it.Set(odf.GetDataPointer()); ++it; } // check if mask image is given if it needs resampling PrepareMaskImage(); // load parameter file - LoadParameters(m_ParameterFile); + LoadParameters(); // prepare parameters float minSpacing; if(m_QBallImage->GetSpacing()[0]GetSpacing()[1] && m_QBallImage->GetSpacing()[0]GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[0]; else if (m_QBallImage->GetSpacing()[1] < m_QBallImage->GetSpacing()[2]) minSpacing = m_QBallImage->GetSpacing()[1]; else minSpacing = m_QBallImage->GetSpacing()[2]; if(m_ParticleLength == 0) m_ParticleLength = 1.5*minSpacing; if(m_ParticleWidth == 0) m_ParticleWidth = 0.5*minSpacing; if(m_ParticleWeight == 0) EstimateParticleWeight(); float alpha = log(m_EndTemperature/m_StartTemperature); m_Steps = m_Iterations/10000; if (m_Steps<10) m_Steps = 10; if (m_Steps>m_Iterations) { MITK_INFO << "GibbsTrackingFilter: not enough iterations!"; m_AbortTracking = true; } if (m_CurvatureThreshold < mitk::eps) m_CurvatureThreshold = 0; unsigned long singleIts = (unsigned long)((1.0*m_Iterations) / (1.0*m_Steps)); // seed random generators Statistics::MersenneTwisterRandomVariateGenerator::Pointer randGen = Statistics::MersenneTwisterRandomVariateGenerator::New(); if (m_RandomSeed>-1) randGen->SetSeed(m_RandomSeed); else randGen->SetSeed(); // load sphere interpolator to evaluate the ODFs SphereInterpolator* interpolator = new SphereInterpolator(m_LutPath); // initialize the actual tracking components (ParticleGrid, Metropolis Hastings Sampler and Energy Computer) ParticleGrid* particleGrid = new ParticleGrid(m_MaskImage, m_ParticleLength); GibbsEnergyComputer* encomp = new GibbsEnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); - -// EnergyComputer* encomp = new EnergyComputer(m_QBallImage, m_MaskImage, particleGrid, interpolator, randGen); encomp->SetParameters(m_ParticleWeight,m_ParticleWidth,m_ConnectionPotential*m_ParticleLength*m_ParticleLength,m_CurvatureThreshold,m_InexBalance,m_ParticlePotential); MetropolisHastingsSampler* sampler = new MetropolisHastingsSampler(particleGrid, encomp, randGen, m_CurvatureThreshold); MITK_INFO << "----------------------------------------"; MITK_INFO << "Iterations: " << m_Iterations; MITK_INFO << "Steps: " << m_Steps; MITK_INFO << "Particle length: " << m_ParticleLength; MITK_INFO << "Particle width: " << m_ParticleWidth; MITK_INFO << "Particle weight: " << m_ParticleWeight; MITK_INFO << "Start temperature: " << m_StartTemperature; MITK_INFO << "End temperature: " << m_EndTemperature; MITK_INFO << "In/Ex balance: " << m_InexBalance; MITK_INFO << "Min. fiber length: " << m_MinFiberLength; MITK_INFO << "Curvature threshold: " << m_CurvatureThreshold; MITK_INFO << "Random seed: " << m_RandomSeed; MITK_INFO << "----------------------------------------"; // main loop + preClock.Stop(); + TimeProbe clock; clock.Start(); m_NumAcceptedFibers = 0; unsigned long counter = 1; for( m_CurrentStep = 1; m_CurrentStep <= m_Steps; m_CurrentStep++ ) { // update temperatur for simulated annealing process float temperature = m_StartTemperature * exp(alpha*(((1.0)*m_CurrentStep)/((1.0)*m_Steps))); sampler->SetTemperature(temperature); for (unsigned long i=0; iMakeProposal(); if (m_BuildFibers || (i==singleIts-1 && m_CurrentStep==m_Steps)) { m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/counter; m_NumParticles = particleGrid->m_NumParticles; m_NumConnections = particleGrid->m_NumConnections; FiberBuilder fiberBuilder(particleGrid, m_MaskImage); m_FiberPolyData = fiberBuilder.iterate(m_MinFiberLength); m_NumAcceptedFibers = m_FiberPolyData->GetNumberOfLines(); m_BuildFibers = false; } counter++; } m_ProposalAcceptance = (float)sampler->GetNumAcceptedProposals()/counter; m_NumParticles = particleGrid->m_NumParticles; m_NumConnections = particleGrid->m_NumConnections; MITK_INFO << "GibbsTrackingFilter: proposal acceptance: " << 100*m_ProposalAcceptance << "%"; MITK_INFO << "GibbsTrackingFilter: particles: " << m_NumParticles; MITK_INFO << "GibbsTrackingFilter: connections: " << m_NumConnections; MITK_INFO << "GibbsTrackingFilter: progress: " << 100*(float)m_CurrentStep/m_Steps << "%"; MITK_INFO << "GibbsTrackingFilter: cell overflows: " << particleGrid->m_NumCellOverflows; MITK_INFO << "----------------------------------------"; if (m_AbortTracking) break; } + clock.Stop(); delete sampler; delete encomp; delete interpolator; delete particleGrid; m_AbortTracking = true; m_BuildFibers = false; - MITK_INFO << "GibbsTrackingFilter: done generate data"; + int h = clock.GetTotal()/3600; + int m = ((int)clock.GetTotal()%3600)/60; + int s = (int)clock.GetTotal()%60; + MITK_INFO << "GibbsTrackingFilter: finished gibbs tracking in " << h << "h, " << m << "m and " << s << "s"; + m = (int)preClock.GetTotal()/60; + s = (int)preClock.GetTotal()%60; + MITK_INFO << "GibbsTrackingFilter: preparation of the data took " << m << "m and " << s << "s"; + MITK_INFO << "GibbsTrackingFilter: " << m_NumAcceptedFibers << " fibers accepted"; + + SaveParameters(); } template< class ItkQBallImageType > void GibbsTrackingFilter< ItkQBallImageType >::PrepareMaskImage() { if(m_MaskImage.IsNull()) { MITK_INFO << "GibbsTrackingFilter: generating default mask image"; m_MaskImage = ItkFloatImageType::New(); m_MaskImage->SetSpacing( m_QBallImage->GetSpacing() ); m_MaskImage->SetOrigin( m_QBallImage->GetOrigin() ); m_MaskImage->SetDirection( m_QBallImage->GetDirection() ); m_MaskImage->SetRegions( m_QBallImage->GetLargestPossibleRegion() ); m_MaskImage->Allocate(); m_MaskImage->FillBuffer(1.0); } else if ( m_MaskImage->GetLargestPossibleRegion().GetSize()[0]!=m_QBallImage->GetLargestPossibleRegion().GetSize()[0] || m_MaskImage->GetLargestPossibleRegion().GetSize()[1]!=m_QBallImage->GetLargestPossibleRegion().GetSize()[1] || m_MaskImage->GetLargestPossibleRegion().GetSize()[2]!=m_QBallImage->GetLargestPossibleRegion().GetSize()[2] || m_MaskImage->GetSpacing()[0]!=m_QBallImage->GetSpacing()[0] || m_MaskImage->GetSpacing()[1]!=m_QBallImage->GetSpacing()[1] || m_MaskImage->GetSpacing()[2]!=m_QBallImage->GetSpacing()[2] ) { MITK_INFO << "GibbsTrackingFilter: resampling mask image"; typedef itk::ResampleImageFilter< ItkFloatImageType, ItkFloatImageType, float > ResamplerType; ResamplerType::Pointer resampler = ResamplerType::New(); resampler->SetOutputSpacing( m_QBallImage->GetSpacing() ); resampler->SetOutputOrigin( m_QBallImage->GetOrigin() ); resampler->SetOutputDirection( m_QBallImage->GetDirection() ); resampler->SetSize( m_QBallImage->GetLargestPossibleRegion().GetSize() ); resampler->SetInput( m_MaskImage ); resampler->SetDefaultPixelValue(1.0); resampler->Update(); m_MaskImage = resampler->GetOutput(); MITK_INFO << "GibbsTrackingFilter: resampling finished"; } } -// load current tracking paramters from xml file (.gtp) +// load tracking paramters from xml file (.gtp) template< class ItkQBallImageType > -bool GibbsTrackingFilter< ItkQBallImageType >::LoadParameters(std::string filename) +bool GibbsTrackingFilter< ItkQBallImageType >::LoadParameters() { m_AbortTracking = true; try { - if( filename.length()==0 ) + if( m_LoadParameterFile.length()==0 ) { m_AbortTracking = false; return true; } - MITK_INFO << "GibbsTrackingFilter: loading parameter file " << filename; + MITK_INFO << "GibbsTrackingFilter: loading parameter file " << m_LoadParameterFile; - TiXmlDocument doc( filename ); + TiXmlDocument doc( m_LoadParameterFile ); doc.LoadFile(); TiXmlHandle hDoc(&doc); TiXmlElement* pElem; TiXmlHandle hRoot(0); pElem = hDoc.FirstChildElement().Element(); hRoot = TiXmlHandle(pElem); pElem = hRoot.FirstChildElement("parameter_set").Element(); QString iterations(pElem->Attribute("iterations")); m_Iterations = iterations.toULong(); QString particleLength(pElem->Attribute("particle_length")); m_ParticleLength = particleLength.toFloat(); QString particleWidth(pElem->Attribute("particle_width")); m_ParticleWidth = particleWidth.toFloat(); QString partWeight(pElem->Attribute("particle_weight")); m_ParticleWeight = partWeight.toFloat(); QString startTemp(pElem->Attribute("temp_start")); m_StartTemperature = startTemp.toFloat(); QString endTemp(pElem->Attribute("temp_end")); m_EndTemperature = endTemp.toFloat(); QString inExBalance(pElem->Attribute("inexbalance")); m_InexBalance = inExBalance.toFloat(); QString fiberLength(pElem->Attribute("fiber_length")); m_MinFiberLength = fiberLength.toFloat(); QString curvThres(pElem->Attribute("curvature_threshold")); m_CurvatureThreshold = cos(curvThres.toFloat()*M_PI/180); m_AbortTracking = false; MITK_INFO << "GibbsTrackingFilter: parameter file loaded successfully"; return true; } catch(...) { MITK_INFO << "GibbsTrackingFilter: could not load parameter file"; return false; } } +// save current tracking paramters to xml file (.gtp) +template< class ItkQBallImageType > +bool GibbsTrackingFilter< ItkQBallImageType >::SaveParameters() +{ + try + { + if( m_SaveParameterFile.length()==0 ) + { + MITK_INFO << "GibbsTrackingFilter: no filename specified to save parameters"; + return true; + } + + MITK_INFO << "GibbsTrackingFilter: saving parameter file " << m_SaveParameterFile; + + TiXmlDocument documentXML; + TiXmlDeclaration* declXML = new TiXmlDeclaration( "1.0", "", "" ); + documentXML.LinkEndChild( declXML ); + + TiXmlElement* mainXML = new TiXmlElement("global_tracking_parameter_file"); + mainXML->SetAttribute("file_version", "0.1"); + documentXML.LinkEndChild(mainXML); + + TiXmlElement* paramXML = new TiXmlElement("parameter_set"); + paramXML->SetAttribute("iterations", QString::number(m_Iterations).toStdString()); + paramXML->SetAttribute("particle_length", QString::number(m_ParticleLength).toStdString()); + paramXML->SetAttribute("particle_width", QString::number(m_ParticleWidth).toStdString()); + paramXML->SetAttribute("particle_weight", QString::number(m_ParticleWeight).toStdString()); + paramXML->SetAttribute("temp_start", QString::number(m_StartTemperature).toStdString()); + paramXML->SetAttribute("temp_end", QString::number(m_EndTemperature).toStdString()); + paramXML->SetAttribute("inexbalance", QString::number(m_InexBalance).toStdString()); + paramXML->SetAttribute("fiber_length", QString::number(m_MinFiberLength).toStdString()); + paramXML->SetAttribute("curvature_threshold", QString::number(m_CurvatureThreshold).toStdString()); + mainXML->LinkEndChild(paramXML); + + QString filename(m_SaveParameterFile.c_str()); + if(!filename.endsWith(".gtp")) + filename += ".gtp"; + documentXML.SaveFile( filename.toStdString() ); + + MITK_INFO << "GibbsTrackingFilter: parameter file saved successfully"; + return true; + } + catch(...) + { + MITK_INFO << "GibbsTrackingFilter: could not save parameter file"; + return false; + } +} + } diff --git a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h index a49d3a3fca..658d1650d5 100644 --- a/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h +++ b/Modules/DiffusionImaging/Tractography/itkGibbsTrackingFilter.h @@ -1,142 +1,145 @@ /*=================================================================== 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 // MITK #include // ITK #include #include #include #include // VTK #include #include #include #include #include namespace itk{ template< class ItkQBallImageType > class GibbsTrackingFilter : public ProcessObject { public: typedef GibbsTrackingFilter Self; typedef ProcessObject Superclass; typedef SmartPointer< Self > Pointer; typedef SmartPointer< const Self > ConstPointer; itkNewMacro(Self) itkTypeMacro( GibbsTrackingFilter, ProcessObject ) typedef Image< DiffusionTensor3D, 3 > ItkTensorImage; typedef typename ItkQBallImageType::Pointer ItkQBallImageTypePointer; typedef Image< float, 3 > ItkFloatImageType; typedef vtkSmartPointer< vtkPolyData > FiberPolyDataType; // parameter setter itkSetMacro( StartTemperature, float ) itkSetMacro( EndTemperature, float ) itkSetMacro( Iterations, unsigned long ) itkSetMacro( ParticleWeight, float ) itkSetMacro( ParticleWidth, float ) itkSetMacro( ParticleLength, float ) itkSetMacro( ConnectionPotential, float ) itkSetMacro( InexBalance, float ) itkSetMacro( ParticlePotential, float ) itkSetMacro( MinFiberLength, int ) itkSetMacro( AbortTracking, bool ) itkSetMacro( CurvatureThreshold, float) itkSetMacro( DuplicateImage, bool ) itkSetMacro( RandomSeed, int ) - itkSetMacro( ParameterFile, std::string ) + itkSetMacro( LoadParameterFile, std::string ) + itkSetMacro( SaveParameterFile, std::string ) itkSetMacro( LutPath, std::string ) // getter itkGetMacro( ParticleWeight, float ) itkGetMacro( ParticleWidth, float ) itkGetMacro( ParticleLength, float ) itkGetMacro( CurrentStep, unsigned long ) itkGetMacro( NumParticles, int ) itkGetMacro( NumConnections, int ) itkGetMacro( NumAcceptedFibers, int ) itkGetMacro( ProposalAcceptance, float ) itkGetMacro( Steps, unsigned int) // input data itkSetMacro(QBallImage, typename ItkQBallImageType::Pointer) itkSetMacro(MaskImage, ItkFloatImageType::Pointer) itkSetMacro(TensorImage, ItkTensorImage::Pointer) void GenerateData(); virtual void Update(){ this->GenerateData(); } FiberPolyDataType GetFiberBundle(); protected: GibbsTrackingFilter(); virtual ~GibbsTrackingFilter(); void EstimateParticleWeight(); void PrepareMaskImage(); - bool LoadParameters(std::string filename); + bool LoadParameters(); + bool SaveParameters(); // Input Images typename ItkQBallImageType::Pointer m_QBallImage; typename ItkFloatImageType::Pointer m_MaskImage; typename ItkTensorImage::Pointer m_TensorImage; // Tracking parameters float m_StartTemperature; // Start temperature float m_EndTemperature; // End temperature unsigned long m_Iterations; // Total number of iterations unsigned long m_CurrentStep; // current tracking step float m_ParticleWeight; // w (unitless) float m_ParticleWidth; // sigma (mm) float m_ParticleLength; // l (mm) float m_ConnectionPotential; // gross L (chemisches potential, default 10) float m_InexBalance; // gewichtung zwischen den lambdas; -5 ... 5 -> nur intern ... nur extern,default 0 float m_ParticlePotential; // default 0.2 int m_MinFiberLength; // discard all fibers shortan than the specified length in mm bool m_AbortTracking; // set flag to abort tracking int m_NumAcceptedFibers; // number of reconstructed fibers generated by the FiberBuilder volatile bool m_BuildFibers; // set flag to generate fibers from particle grid unsigned int m_Steps; // number of temperature decrease steps float m_ProposalAcceptance; // proposal acceptance rate (0-1) float m_CurvatureThreshold; // curvature threshold in radians (1 -> no curvature is accepted, -1 all curvature angles are accepted) bool m_DuplicateImage; // generates a working copy of the qball image so that the original image won't be changed by the mean subtraction int m_NumParticles; // current number of particles in grid int m_NumConnections; // current number of connections between particles in grid int m_RandomSeed; // seed value for random generator (-1 for standard seeding) - std::string m_ParameterFile; // filename of parameter file + std::string m_LoadParameterFile; // filename of parameter file (reader) + std::string m_SaveParameterFile; // filename of parameter file (writer) std::string m_LutPath; // path to lookuptables used by the sphere interpolator FiberPolyDataType m_FiberPolyData; // container for reconstructed fibers }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkGibbsTrackingFilter.cpp" #endif #endif