diff --git a/Modules/PhotoacousticsLib/include/mitkPAVessel.h b/Modules/PhotoacousticsLib/include/mitkPAVessel.h index 8670e8c684..25a223df17 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVessel.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVessel.h @@ -1,117 +1,117 @@ /*=================================================================== 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 "mitkPAVesselMeanderStrategy.h" #include "mitkPAInSilicoTissueVolume.h" #include "mitkPAVector.h" #include "mitkPAVesselProperties.h" #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT Vessel : public itk::LightObject { public: mitkClassMacroItkParent(Vessel, itk::LightObject) mitkNewMacro1Param(Self, VesselProperties::Pointer) /** * Callback function definition of a VesselMeanderStrategy */ typedef void (VesselMeanderStrategy::*CalculateNewVesselPositionCallback) (Vector::Pointer, Vector::Pointer, double, std::mt19937*); /** * @brief ExpandVessel makes this Vessel expand one step in its current direction. * After expanding, the vessel will draw itself into the given InSilicoTissueVolume. * * @param volume volume for the vessel to draw itself in * @param calculateNewPosition a callback function of the VesselMeanderStrategy class. * It is used to calculate the final position after taking the step. * @param bendingFactor a metric of how much the Vessel should bend. If set to 0 the vessel will go in a straight line. */ void ExpandVessel(mitk::pa::InSilicoTissueVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng); /** * @brief CanBifurcate * @return true if the Vessel 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. */ Vessel::Pointer Bifurcate(std::mt19937* rng); /** * @brief IsFinished * @return true if the vessel cannot expand any further */ bool IsFinished(); itkGetConstMacro(VesselProperties, VesselProperties::Pointer); protected: Vessel(VesselProperties::Pointer parameters); virtual ~Vessel(); private: - const double MINIMUM_VESSEL_RADIUS = 1; + const double MINIMUM_VESSEL_RADIUS = 0.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(Vector::Pointer toPosition, mitk::pa::InSilicoTissueVolume::Pointer volume); VesselProperties::Pointer m_VesselProperties; VesselMeanderStrategy::Pointer m_VesselMeanderStrategy; bool m_Finished; double m_WalkedDistance; std::uniform_real_distribution<> m_RangeDistribution; std::uniform_real_distribution<> m_SignDistribution; std::uniform_real_distribution<> m_RadiusRangeDistribution; int GetSign(std::mt19937* rng); }; /** * @brief Equal A function comparing two vessels for beeing equal * * @param rightHandSide A vessel to be compared * @param leftHandSide A vessel to be compared * @param eps tolarence for comparison. You can use mitk::eps in most cases. * @param verbose flag indicating if the user wants detailed console output or not. * @return true, if all subsequent comparisons are true, false otherwise */ MITKPHOTOACOUSTICSLIB_EXPORT bool Equal(const Vessel::Pointer leftHandSide, const Vessel::Pointer rightHandSide, double eps, bool verbose); } } #endif // MITKVESSEL_H diff --git a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp index c33c71618c..f785185151 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp @@ -1,169 +1,206 @@ /*=================================================================== 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 "mitkPAVessel.h" #include #define _USE_MATH_DEFINES #include #include mitk::pa::Vessel::Vessel(VesselProperties::Pointer initialProperties) : m_RangeDistribution(M_PI / 16, M_PI / 8), m_SignDistribution(-1, 1) { m_Finished = false; //Copy this so it may be reused for other vessels. m_VesselProperties = VesselProperties::New(initialProperties); m_RadiusRangeDistribution = std::uniform_real_distribution<>(NEW_RADIUS_MINIMUM_RELATIVE_SIZE, NEW_RADIUS_MAXIMUM_RELATIVE_SIZE); m_VesselMeanderStrategy = VesselMeanderStrategy::New(); m_WalkedDistance = 0; } mitk::pa::Vessel::~Vessel() { m_VesselProperties = nullptr; m_VesselMeanderStrategy = nullptr; } void mitk::pa::Vessel::ExpandVessel(InSilicoTissueVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng) { Vector::Pointer oldPosition = m_VesselProperties->GetPositionVector()->Clone(); (m_VesselMeanderStrategy->*calculateNewPosition)(m_VesselProperties->GetPositionVector(), m_VesselProperties->GetDirectionVector(), bendingFactor, rng); DrawVesselInVolume(oldPosition, volume); } bool mitk::pa::Vessel::CanBifurcate() { return m_VesselProperties->GetBifurcationFrequency() < m_WalkedDistance; } int mitk::pa::Vessel::GetSign(std::mt19937 *rng) { if (m_SignDistribution(*rng) < 0) return -1; return 1; } mitk::pa::Vessel::Pointer mitk::pa::Vessel::Bifurcate(std::mt19937* rng) { VesselProperties::Pointer vesselParams = VesselProperties::New(m_VesselProperties); double thetaChange = m_RangeDistribution(*rng) * GetSign(rng); double phiChange = m_RangeDistribution(*rng) * GetSign(rng); vesselParams->GetDirectionVector()->Rotate(thetaChange, phiChange); m_VesselProperties->GetDirectionVector()->Rotate(-thetaChange, -phiChange); double newRadius = m_RadiusRangeDistribution(*rng)*m_VesselProperties->GetRadiusInVoxel(); vesselParams->SetRadiusInVoxel(newRadius); m_VesselProperties->SetRadiusInVoxel( sqrt(m_VesselProperties->GetRadiusInVoxel()*m_VesselProperties->GetRadiusInVoxel() - newRadius*newRadius)); m_WalkedDistance = 0; return Vessel::New(vesselParams); } void mitk::pa::Vessel::DrawVesselInVolume(Vector::Pointer fromPosition, InSilicoTissueVolume::Pointer volume) { Vector::Pointer diffVector = Vector::New(); Vector::Pointer toPosition = m_VesselProperties->GetPositionVector(); diffVector->SetElement(0, fromPosition->GetElement(0) - toPosition->GetElement(0)); diffVector->SetElement(1, fromPosition->GetElement(1) - toPosition->GetElement(1)); diffVector->SetElement(2, fromPosition->GetElement(2) - toPosition->GetElement(2)); //1/SCALING_FACTOR steps along the direction vector are taken and drawn into the image. Vector::Pointer stepSize = Vector::New(); stepSize->SetValue(m_VesselProperties->GetDirectionVector()); 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); + double xPos = fromPosition->GetElement(0); + double yPos = fromPosition->GetElement(1); + double zPos = fromPosition->GetElement(2); if (!volume->IsInsideVolume(xPos, yPos, zPos)) { m_VesselProperties->SetRadiusInVoxel(0); break; } double radius = m_VesselProperties->GetRadiusInVoxel(); + double ceiledRadius = ceil(radius); - for (int x = xPos - radius; x <= xPos + radius; x++) - for (int y = yPos - radius; y <= yPos + radius; y++) - for (int z = zPos - radius; z <= zPos + radius; z++) + for (int x = xPos - ceiledRadius; x <= xPos + ceiledRadius; x += 1) + for (int y = yPos - ceiledRadius; y <= yPos + ceiledRadius; y += 1) + for (int z = zPos - ceiledRadius; z <= zPos + ceiledRadius; z += 1) { - if (radius*radius >= (x - xPos)*(x - xPos) + (y - yPos)*(y - yPos) + (z - zPos)*(z - zPos)) + if (!volume->IsInsideVolume(x, y, z)) { - volume->SetVolumeValues(x, y, z, m_VesselProperties->GetAbsorptionCoefficient(), + continue; + } + double xDiff = x - xPos; + double yDiff = y - yPos; + double zDiff = z - zPos; + double vectorLengthDiff = radius*radius - (xDiff*xDiff + yDiff*yDiff + zDiff*zDiff); + + if (vectorLengthDiff > 0) + { + volume->SetVolumeValues(x, y, z, + m_VesselProperties->GetAbsorptionCoefficient(), m_VesselProperties->GetScatteringCoefficient(), m_VesselProperties->GetAnisotopyCoefficient(), mitk::pa::InSilicoTissueVolume::SegmentationType::VESSEL); } + else + { + double backgroundFraction = abs(sqrt((xDiff*xDiff + yDiff*yDiff + zDiff*zDiff)) - radius); + double vesselFraction = 1.0 - backgroundFraction; + auto type = mitk::pa::InSilicoTissueVolume::SegmentationType::BACKGROUND; + if (vesselFraction >= 0.5) + { + type = mitk::pa::InSilicoTissueVolume::SegmentationType::VESSEL; + } + double absorption = backgroundFraction * volume->GetAbsorptionVolume()->GetData( + x, y, z) + + vesselFraction * m_VesselProperties->GetAbsorptionCoefficient(); + double scattering = backgroundFraction * volume->GetScatteringVolume()->GetData( + x, y, z) + + vesselFraction * m_VesselProperties->GetScatteringCoefficient(); + double anisotropy = backgroundFraction * volume->GetAnisotropyVolume()->GetData( + x, y, z) + + vesselFraction * m_VesselProperties->GetAnisotopyCoefficient(); + + volume->SetVolumeValues(x, y, z, + absorption, + scattering, + anisotropy, + type + ); + } } diffVector->SetElement(0, fromPosition->GetElement(0) - toPosition->GetElement(0)); diffVector->SetElement(1, fromPosition->GetElement(1) - toPosition->GetElement(1)); diffVector->SetElement(2, fromPosition->GetElement(2) - toPosition->GetElement(2)); } } bool mitk::pa::Vessel::IsFinished() { return m_VesselProperties->GetRadiusInVoxel() < MINIMUM_VESSEL_RADIUS; } bool mitk::pa::Equal(const Vessel::Pointer leftHandSide, const Vessel::Pointer rightHandSide, double eps, bool verbose) { MITK_INFO(verbose) << "=== mitk::pa::Vessel Equal ==="; if (rightHandSide.IsNull() || leftHandSide.IsNull()) { MITK_INFO(verbose) << "Cannot compare nullpointers"; return false; } if (leftHandSide->IsFinished() != rightHandSide->IsFinished()) { MITK_INFO(verbose) << "Not same finished state."; return false; } if (leftHandSide->CanBifurcate() != rightHandSide->CanBifurcate()) { MITK_INFO(verbose) << "Not same bifurcation state."; return false; } if (!Equal(leftHandSide->GetVesselProperties(), rightHandSide->GetVesselProperties(), eps, verbose)) { MITK_INFO(verbose) << "Vesselproperties not equal"; return false; } return true; }