diff --git a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp index 7caf51845c..cfe28c043c 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp @@ -1,217 +1,253 @@ /*=================================================================== 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)); + Volume::Pointer originalVolumeAbsorption = volume->GetAbsorptionVolume()->DeepCopy(); + Volume::Pointer originalVolumeScattering = volume->GetScatteringVolume()->DeepCopy(); + Volume::Pointer originalVolumeAnisotropy = volume->GetAnisotropyVolume()->DeepCopy(); + //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)); 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); + double radiusSquared = radius*radius; + double PI = 3.141; + unsigned int totalPoints = 8000; 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 (!volume->IsInsideVolume(x, y, z)) { continue; } double xDiff = x - xPos; double yDiff = y - yPos; double zDiff = z - zPos; - double vectorLengthDiff = radius*radius - (xDiff*xDiff + yDiff*yDiff + zDiff*zDiff); + double vectorLengthDiff = radius - sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff); + + if (vectorLengthDiff < -2) + { + continue; + } - if (vectorLengthDiff > 0) + if (vectorLengthDiff > 1) { volume->SetVolumeValues(x, y, z, m_VesselProperties->GetAbsorptionCoefficient(), m_VesselProperties->GetScatteringCoefficient(), m_VesselProperties->GetAnisotopyCoefficient(), mitk::pa::InSilicoTissueVolume::SegmentationType::VESSEL); } else { - if (m_VesselProperties->GetDoPartialVolume()) + if (m_VesselProperties->GetDoPartialVolume() && volume->GetSegmentationVolume()->GetData(x, y, z) != InSilicoTissueVolume::SegmentationType::VESSEL) { - double backgroundFraction = abs(sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff) - radius); - if (backgroundFraction <= 1 && backgroundFraction >= 0) + unsigned int pointsWithin = 0; + + double phi = atan2(y, x); + double theta = atan2(sqrt(x*x + y*y), z); + double signx = (phi < PI || phi >(3 * PI / 2)) ? -1 : 1; + double signy = (phi < PI) ? -1 : 1; + double signz = (theta < PI / 2) ? -1 : 1; + + for (double checkx = 0; checkx < 1; checkx += 0.05) + for (double checky = 0; checky < 1; checky += 0.05) + for (double checkz = 0; checkz < 1; checkz += 0.05) + { + double checkxSign = checkx*signx; + double checkySign = checky*signy; + double checkzSign = checkz*signz; + if (((xDiff - checkxSign) + (xDiff - checkxSign) + + (yDiff - checkySign) + (yDiff - checkySign) + + (zDiff - checkzSign) + (zDiff - checkzSign)) < radiusSquared) + { + pointsWithin += 1; + } + } + + double vesselFraction = pointsWithin / totalPoints; + + if (vesselFraction <= 1 && vesselFraction >= 0) { - double vesselFraction = 1.0 - backgroundFraction; + double backgroundFraction = 1.0 - vesselFraction; - double absorption = backgroundFraction * volume->GetAbsorptionVolume()->GetData( + double absorption = backgroundFraction * originalVolumeAbsorption->GetData( x, y, z) + vesselFraction * m_VesselProperties->GetAbsorptionCoefficient(); - double scattering = backgroundFraction * volume->GetScatteringVolume()->GetData( + double scattering = backgroundFraction * originalVolumeScattering->GetData( x, y, z) + vesselFraction * m_VesselProperties->GetScatteringCoefficient(); - double anisotropy = backgroundFraction * volume->GetAnisotropyVolume()->GetData( + double anisotropy = backgroundFraction * originalVolumeAnisotropy->GetData( x, y, z) + vesselFraction * m_VesselProperties->GetAnisotopyCoefficient(); if (vesselFraction >= 0.5) { volume->SetVolumeValues(x, y, z, absorption, scattering, anisotropy, InSilicoTissueVolume::SegmentationType::VESSEL); } else { volume->SetVolumeValues(x, y, z, absorption, scattering, anisotropy); } } } } } 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; } diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp index b319b94bb3..0306de205f 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp @@ -1,169 +1,219 @@ /*=================================================================== 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 "mitkPAInSilicoTissueVolume.h" #include "mitkPAVector.h" #include "mitkPAVessel.h" #include "mitkIOUtil.h" class mitkPhotoacousticVesselTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticVesselTestSuite); MITK_TEST(testEmptyInitializationProperties); MITK_TEST(testWalkInStraightLine); MITK_TEST(testBifurcate); + MITK_TEST(testPartialVolume); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::Vessel::Pointer m_TestVessel; mitk::pa::Vessel::CalculateNewVesselPositionCallback m_StraightLine; mitk::pa::Vessel::CalculateNewVesselPositionCallback m_Diverging; mitk::pa::InSilicoTissueVolume::Pointer m_TestInSilicoVolume; mitk::pa::TissueGeneratorParameters::Pointer m_TestVolumeParameters; public: void setUp() { auto params = mitk::pa::VesselProperties::New(); m_TestVessel = mitk::pa::Vessel::New(params); m_StraightLine = &mitk::pa::VesselMeanderStrategy::CalculateNewPositionInStraightLine; m_Diverging = &mitk::pa::VesselMeanderStrategy::CalculateRandomlyDivergingPosition; m_TestVolumeParameters = createTestVolumeParameters(); m_TestInSilicoVolume = mitk::pa::InSilicoTissueVolume::New(m_TestVolumeParameters); } mitk::pa::TissueGeneratorParameters::Pointer createTestVolumeParameters() { auto returnParameters = mitk::pa::TissueGeneratorParameters::New(); returnParameters->SetXDim(10); returnParameters->SetYDim(10); returnParameters->SetZDim(10); returnParameters->SetBackgroundAbsorption(0); returnParameters->SetBackgroundScattering(0); returnParameters->SetBackgroundAnisotropy(0); return returnParameters; } void testEmptyInitializationProperties() { CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == true); } + void testPartialVolume() + { + auto testPosition = mitk::pa::Vector::New(); + testPosition->SetElement(0, 0); + testPosition->SetElement(1, 4); + testPosition->SetElement(2, 4); + auto testDirection = mitk::pa::Vector::New(); + testDirection->SetElement(0, 1); + testDirection->SetElement(1, 0); + testDirection->SetElement(2, 0); + auto params = mitk::pa::VesselProperties::New(); + params->SetDoPartialVolume(true); + params->SetRadiusInVoxel(1); + params->SetBifurcationFrequency(100); + params->SetAbsorptionCoefficient(10); + params->SetScatteringCoefficient(10); + params->SetAnisotopyCoefficient(10); + params->SetPositionVector(testPosition); + params->SetDirectionVector(testDirection); + m_TestVessel = mitk::pa::Vessel::New(params); + + CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); + CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); + + CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); + + m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); + + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 3)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 3) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 4)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 5)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 5) - 5.85786438) <= 1e-6); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 3)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 3) - 5.85786438) <= 1e-6); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 3)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 3) - 5.85786438) <= 1e-6); + CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 5)), + abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 5) - 5.85786438) <= 1e-6); + } + void testWalkInStraightLine() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); + params->SetDoPartialVolume(false); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(100); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 6, 4)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5) - 10) <= mitk::eps); + CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 6)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 5, 4) - 10) <= mitk::eps); + CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 5, 4)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 6, 4)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); - CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 5) - 10) <= mitk::eps); + CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 5)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 6)) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(3, 4, 4)) <= mitk::eps); } void testBifurcate() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(1); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == true); std::mt19937 rng; auto bifurcationVessel = m_TestVessel->Bifurcate(&rng); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(bifurcationVessel->CanBifurcate() == false); } void tearDown() { } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVessel)