diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp index 432423ed99..d5ac1021de 100644 --- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp +++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp @@ -1,163 +1,167 @@ /*=================================================================== 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 "mitkPhotoacousticVessel.h" #include #include #include mitk::PhotoacousticVessel::PhotoacousticVessel(): - m_Range(M_PI/16, M_PI/8), - m_Sign(-1, 1) + m_RangeDistribution(M_PI/16, M_PI/8), + m_SignDistribution(-1, 1) { srand(std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count()); m_Finished = false; m_Position = nullptr; m_Direction = nullptr; m_Radius = -1; m_AbsorptionCoefficient = -1; m_ScatteringCoefficient = -1; m_AnisotropyCoefficient = -1; m_VesselMeanderStrategy = nullptr; m_WalkedDistance = 0; m_State = UNINITIALIZED; } mitk::PhotoacousticVessel::~PhotoacousticVessel() { m_Position = nullptr; m_Direction = nullptr; m_VesselMeanderStrategy = nullptr; m_State = UNINITIALIZED; } void mitk::PhotoacousticVessel::Initialize(mitk::PhotoacousticSmartVector::Pointer initialPosition, mitk::PhotoacousticSmartVector::Pointer initialDirection, double initialRadius, double absorptionCoefficient, double scatteringCoefficient, double anisotophyCoefficient, double bifurcationFrequency, std::mt19937* rng) { m_Radius = initialRadius; m_AbsorptionCoefficient = absorptionCoefficient; m_ScatteringCoefficient = scatteringCoefficient; m_AnisotropyCoefficient = anisotophyCoefficient; m_Finished = false; m_Position = initialPosition; m_Direction = initialDirection; m_VesselMeanderStrategy = mitk::PhotoacousticVesselMeanderStrategy::New(); m_State = INITIALIZED; - std::normal_distribution<> factor(1, 0.2); + std::normal_distribution<> factor(1, VESSEL_BIFURCATION_VARIATION_SIGMA); m_BifurcationFrequency = bifurcationFrequency * factor(*rng); - MITK_INFO << "Bif. frequency: " << m_BifurcationFrequency; + m_RadiusRangeDistribution = std::uniform_real_distribution<>(m_Radius * NEW_RADIUS_MINIMUM_RELATIVE_SIZE, + m_Radius * NEW_RADIUS_MAXIMUM_RELATIVE_SIZE); } void mitk::PhotoacousticVessel::ExpandVessel(mitk::PhotoacousticVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng) { AssertState(INITIALIZED); mitk::PhotoacousticSmartVector::Pointer oldPosition = mitk::PhotoacousticSmartVector::New(); oldPosition->SetValue(m_Position); (m_VesselMeanderStrategy->*calculateNewPosition)(m_Position, m_Direction, bendingFactor, rng); DrawVesselInVolume(oldPosition, volume); } bool mitk::PhotoacousticVessel::CanBifurcate() { return m_BifurcationFrequency < m_WalkedDistance; } +int mitk::PhotoacousticVessel::GetSign(std::mt19937 *rng) +{ + if(m_SignDistribution(*rng) < 0) + return -1; + + return 1; +} + mitk::PhotoacousticVessel::Pointer mitk::PhotoacousticVessel::Bifurcate(std::mt19937* rng) { AssertState(INITIALIZED); mitk::PhotoacousticVessel::Pointer newVessel = mitk::PhotoacousticVessel::New(); mitk::PhotoacousticSmartVector::Pointer newDirection = m_Direction->Clone(); - double thetaChange = m_Range(*rng) * (m_Sign(*rng) < 0 ? -1 : 1); - double phiChange = m_Range(*rng) * (m_Sign(*rng) < 0 ? -1 : 1); - - MITK_DEBUG << "Theta change: " << thetaChange; - MITK_DEBUG << "Phi change: " << phiChange; - MITK_DEBUG << "(m_Sign(*rng) = " << (m_Sign(*rng)); + double thetaChange = m_RangeDistribution(*rng) * GetSign(rng); + double phiChange = m_RangeDistribution(*rng) * GetSign(rng); newDirection->Rotate(thetaChange, phiChange); m_Direction->Rotate(-thetaChange, -phiChange); - std::uniform_real_distribution<> rangeRadius(m_Radius*0.6, m_Radius*0.8); - - double newRadius = rangeRadius(*rng); + double newRadius = m_RadiusRangeDistribution(*rng); m_Radius = sqrt(m_Radius*m_Radius - newRadius*newRadius); newVessel->Initialize(m_Position->Clone(), newDirection, newRadius, m_AbsorptionCoefficient, m_ScatteringCoefficient, m_AnisotropyCoefficient, m_BifurcationFrequency, rng); m_WalkedDistance = 0; return newVessel; } void mitk::PhotoacousticVessel::DrawVesselInVolume(mitk::PhotoacousticSmartVector::Pointer fromPosition, mitk::PhotoacousticVolume::Pointer volume) { AssertState(INITIALIZED); mitk::PhotoacousticSmartVector::Pointer diffVector = mitk::PhotoacousticSmartVector::New(); diffVector->SetElement(0, fromPosition->GetElement(0)-m_Position->GetElement(0)); diffVector->SetElement(1, fromPosition->GetElement(1)-m_Position->GetElement(1)); diffVector->SetElement(2, fromPosition->GetElement(2)-m_Position->GetElement(2)); + //1/SCALING_FACTOR steps along the direction vector are taken and drawn into the image. mitk::PhotoacousticSmartVector::Pointer stepSize = mitk::PhotoacousticSmartVector::New(); stepSize->SetValue(m_Direction); stepSize->Scale(SCALING_FACTOR); while (diffVector->GetNorm() >= SCALING_FACTOR) { m_WalkedDistance += stepSize->GetNorm(); fromPosition->SetElement(0, fromPosition->GetElement(0)+stepSize->GetElement(0)); fromPosition->SetElement(1, fromPosition->GetElement(1)+stepSize->GetElement(1)); fromPosition->SetElement(2, fromPosition->GetElement(2)+stepSize->GetElement(2)); int xPos = fromPosition->GetElement(0); int yPos = fromPosition->GetElement(1); int zPos = fromPosition->GetElement(2); for(int x = xPos - m_Radius; x <= xPos + m_Radius; x++) for(int y = yPos - m_Radius; y <= yPos + m_Radius; y++) for(int z = zPos - m_Radius; z <= zPos + m_Radius; z++) { if( m_Radius*m_Radius >= (x-xPos)*(x-xPos)+(y-yPos)*(y-yPos)+(z-zPos)*(z-zPos)) { volume->SetVolumeValues(x, y, z, m_AbsorptionCoefficient, m_ScatteringCoefficient, m_AnisotropyCoefficient, mitk::PhotoacousticVolume::SegmentationType::VESSEL); } } diffVector->SetElement(0, fromPosition->GetElement(0)-m_Position->GetElement(0)); diffVector->SetElement(1, fromPosition->GetElement(1)-m_Position->GetElement(1)); diffVector->SetElement(2, fromPosition->GetElement(2)-m_Position->GetElement(2)); if(!volume->IsInsideVolume(xPos, yPos, zPos)) { m_Radius = 0; break; } } } bool mitk::PhotoacousticVessel::IsFinished() { - return m_Radius < MIN_RADIUS; + return m_Radius < MINIMUM_VESSEL_RADIUS; } diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h index 0e6364acfd..10d5a990a4 100644 --- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h +++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h @@ -1,119 +1,125 @@ /*=================================================================== 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 MITKVESSEL_H #define MITKVESSEL_H #include "mitkVector.h" #include "mitkPhotoacousticVesselMeanderStrategy.h" #include "mitkPhotoacousticVolume.h" #include "mitkPhotoacousticSmartVector.h" #include "mitkPhotoacousticStatefulObject.h" #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticVessel : public itk::LightObject, public mitk::PhotoacousticStatefulObject { public: mitkClassMacroItkParent(mitk::PhotoacousticVessel, itk::LightObject) itkFactorylessNewMacro(Self) /** * Callback function definition of a PhotoacousticVesselMeanderStrategy */ typedef void (mitk::PhotoacousticVesselMeanderStrategy::*CalculateNewVesselPositionCallback) (mitk::PhotoacousticSmartVector::Pointer, mitk::PhotoacousticSmartVector::Pointer, double, std::mt19937*); /** * @brief Initialize initializes this PhotoacousticVessel with the given parameters. * * @param initialPosition * @param initialDirection * @param initialRadius * @param absorptionCoefficient * @param scatteringCoefficient * @param anisotropyCoefficient * @param bifurcationFrequency * @param rng */ void Initialize(mitk::PhotoacousticSmartVector::Pointer initialPosition, mitk::PhotoacousticSmartVector::Pointer initialDirection, double initialRadius, double absorptionCoefficient, double scatteringCoefficient, double anisotropyCoefficient, double bifurcationFrequency, std::mt19937* rng); /** * @brief ExpandVessel makes this PhotoacousticVessel expand one step in its current direction. * After expanding, the vessel will draw itself into the given PhotoacousticVolume. * * @param volume volume for the vessel to draw itself in * @param calculateNewPosition a callback function of the PhotoacousticVesselMeanderStrategy class. * It is used to calculate the final position after taking the step. * @param bendingFactor a metric of how much the PhotoacousticVessel should bend. If set to 0 the vessel will go in a straight line. */ void ExpandVessel(mitk::PhotoacousticVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng); /** * @brief CanBifurcate * @return true if the PhotoacousticVessel is ready to Bifurcate() */ bool CanBifurcate(); /** * @brief Bifurcate bifurcates this vessel into two new ones. Makes sure that the volume of the vessels stays the same. * * @return a new vessel split up from the current one. */ mitk::PhotoacousticVessel::Pointer Bifurcate(std::mt19937* rng); /** * @brief IsFinished * @return true if the vessel cannot expand any further */ bool IsFinished(); protected: PhotoacousticVessel(); virtual ~PhotoacousticVessel(); private: const double MIN_RADIUS = 1; const double SCALING_FACTOR = 0.33; + const double VESSEL_BIFURCATION_VARIATION_SIGMA = 0.2; + const double NEW_RADIUS_MINIMUM_RELATIVE_SIZE = 0.6; + const double NEW_RADIUS_MAXIMUM_RELATIVE_SIZE = 0.8; void DrawVesselInVolume(mitk::PhotoacousticSmartVector::Pointer toPosition, mitk::PhotoacousticVolume::Pointer volume); mitk::PhotoacousticSmartVector::Pointer m_Position; mitk::PhotoacousticSmartVector::Pointer m_Direction; double m_Radius; double m_AbsorptionCoefficient; double m_ScatteringCoefficient; double m_AnisotropyCoefficient; double m_WalkedDistance; double m_BifurcationFrequency; mitk::PhotoacousticVesselMeanderStrategy::Pointer m_VesselMeanderStrategy; bool m_Finished; - std::uniform_real_distribution<> m_Range; - std::uniform_real_distribution<> m_Sign; + std::uniform_real_distribution<> m_RangeDistribution; + std::uniform_real_distribution<> m_SignDistribution; + std::uniform_real_distribution<> m_RadiusRangeDistribution; + + int GetSign(std::mt19937* rng); }; } #endif // MITKVESSEL_H