diff --git a/Modules/PhotoacousticsLib/files.cmake b/Modules/PhotoacousticsLib/files.cmake index 8b14417421..28e6b805ad 100644 --- a/Modules/PhotoacousticsLib/files.cmake +++ b/Modules/PhotoacousticsLib/files.cmake @@ -1,52 +1,54 @@ SET(H_FILES include/mitkPAPropertyCalculator.h include/mitkPAVector.h include/mitkPATissueGeneratorParameters.h include/mitkPAInSilicoTissueVolume.h include/mitkPATissueGenerator.h include/mitkPAVesselTree.h include/mitkPAVessel.h + include/mitkPAVesselDrawer.h include/mitkPAVesselMeanderStrategy.h include/mitkPANoiseGenerator.h include/mitkPAVolume.h include/mitkPAComposedVolume.h include/mitkPASlicedVolumeGenerator.h include/mitkPAProbe.h include/mitkPALightSource.h include/mitkPAIOUtil.h include/mitkPAMonteCarloThreadHandler.h include/mitkPASimulationBatchGenerator.h include/mitkPAFluenceYOffsetPair.h include/mitkPAVolumeManipulator.h include/mitkPAVesselProperties.h include/mitkPASimulationBatchGeneratorParameters.h include/mitkPAExceptions.h ) set(CPP_FILES Domain/Vessel/mitkPAVesselTree.cpp Domain/Vessel/mitkPAVessel.cpp Domain/Vessel/mitkPAVesselMeanderStrategy.cpp Domain/Vessel/mitkPAVesselProperties.cpp Domain/Volume/mitkPAInSilicoTissueVolume.cpp Domain/Volume/mitkPAVolume.cpp Domain/Volume/mitkPAComposedVolume.cpp Domain/Volume/mitkPAFluenceYOffsetPair.cpp Generator/mitkPATissueGenerator.cpp Generator/mitkPANoiseGenerator.cpp Generator/mitkPASlicedVolumeGenerator.cpp Generator/mitkPASimulationBatchGenerator.cpp Generator/mitkPASimulationBatchGeneratorParameters.cpp IO/mitkPAIOUtil.cpp Utils/mitkPAPropertyCalculator.cpp Utils/mitkPAVector.cpp Utils/mitkPATissueGeneratorParameters.cpp Utils/mitkPAVolumeManipulator.cpp Utils/ProbeDesign/mitkPAProbe.cpp Utils/ProbeDesign/mitkPALightSource.cpp Utils/Thread/mitkPAMonteCarloThreadHandler.cpp + Utils/mitkPAVesselDrawer.cpp ) set(RESOURCE_FILES spectralLIB.dat ) diff --git a/Modules/PhotoacousticsLib/include/mitkPAVessel.h b/Modules/PhotoacousticsLib/include/mitkPAVessel.h index 25a223df17..2e9f35ee3b 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVessel.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVessel.h @@ -1,117 +1,119 @@ /*=================================================================== 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 "mitkPAVesselDrawer.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 = 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); + + VesselProperties::Pointer m_VesselProperties; + + VesselDrawer::Pointer m_VesselDrawer; }; /** * @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/include/mitkPAVesselDrawer.h b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h new file mode 100644 index 0000000000..3e0f225ec6 --- /dev/null +++ b/Modules/PhotoacousticsLib/include/mitkPAVesselDrawer.h @@ -0,0 +1,82 @@ +/*=================================================================== + +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 MITKVESSELDRAWER_H +#define MITKVESSELDRAWER_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 VesselDrawer : public itk::LightObject + { + public: + mitkClassMacroItkParent(VesselDrawer, itk::LightObject); + itkFactorylessNewMacro(Self); + + void DrawVesselInVolume(Vector::Pointer toPosition, + VesselProperties::Pointer properties, + InSilicoTissueVolume::Pointer volume); + + itkGetMacro(VolumeState, Volume::Pointer); + + protected: + VesselDrawer(); + virtual ~VesselDrawer(); + + private: + + void InitializeVolumeState(InSilicoTissueVolume::Pointer volume); + + enum VolumeState + { + UNTOUCHED, PARTIAL, FULL + }; + + Volume::Pointer m_VolumeState; + Volume::Pointer m_OriginalVolumeAbsorption; + Volume::Pointer m_OriginalVolumeScattering; + Volume::Pointer m_OriginalVolumeAnisotropy; + + void FillInVolumeValues(mitk::pa::InSilicoTissueVolume::Pointer volume, + mitk::pa::VesselProperties::Pointer properties, + double vesselFraction, + int x, int y, int z) + }; + + /** + * @brief Equal A function comparing two vessel drawer instances for beeing equal + * + * @param rightHandSide A vesseldrawer to be compared + * @param leftHandSide A vesseldrawer 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 VesselDrawer::Pointer leftHandSide, const VesselDrawer::Pointer rightHandSide, double eps, bool verbose); + } +} +#endif // MITKVESSELDRAWER_H diff --git a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp index cfe28c043c..09cc8dc1f3 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp @@ -1,253 +1,118 @@ /*=================================================================== 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; + m_VesselDrawer = VesselDrawer::New(); } 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); + m_VesselDrawer->DrawVesselInVolume(oldPosition, m_VesselProperties, volume); + m_WalkedDistance += m_VesselProperties->GetDirectionVector()->GetNorm(); } 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 - sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff); - - if (vectorLengthDiff < -2) - { - continue; - } - - 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() && volume->GetSegmentationVolume()->GetData(x, y, z) != InSilicoTissueVolume::SegmentationType::VESSEL) - { - 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 backgroundFraction = 1.0 - vesselFraction; - - double absorption = backgroundFraction * originalVolumeAbsorption->GetData( - x, y, z) - + vesselFraction * m_VesselProperties->GetAbsorptionCoefficient(); - double scattering = backgroundFraction * originalVolumeScattering->GetData( - x, y, z) - + vesselFraction * m_VesselProperties->GetScatteringCoefficient(); - 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/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp index 699235bdac..679f0bda4a 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVesselMeanderStrategy.cpp @@ -1,53 +1,59 @@ /*=================================================================== 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 "mitkPAVesselMeanderStrategy.h" mitk::pa::VesselMeanderStrategy::VesselMeanderStrategy() { } mitk::pa::VesselMeanderStrategy::~VesselMeanderStrategy() { } void mitk::pa::VesselMeanderStrategy::CalculateNewPositionInStraightLine( Vector::Pointer position, Vector::Pointer direction, double /*bendingFactor*/, std::mt19937* rng) { if (direction->GetNorm() <= mitk::eps) { direction->Randomize(rng); } + direction->Normalize(); + position->SetElement(0, position->GetElement(0) + direction->GetElement(0)); position->SetElement(1, position->GetElement(1) + direction->GetElement(1)); position->SetElement(2, position->GetElement(2) + direction->GetElement(2)); + + position->Normalize(); } void mitk::pa::VesselMeanderStrategy::CalculateRandomlyDivergingPosition( Vector::Pointer position, Vector::Pointer direction, double bendingFactor, std::mt19937* rng) { if (direction->GetNorm() <= mitk::eps) { direction->Randomize(rng); } direction->RandomizeByPercentage(RANDOMIZATION_PERCENTAGE, bendingFactor, rng); + direction->Normalize(); + position->SetElement(0, position->GetElement(0) + direction->GetElement(0)); position->SetElement(1, position->GetElement(1) + direction->GetElement(1)); position->SetElement(2, position->GetElement(2) + direction->GetElement(2)); } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp new file mode 100644 index 0000000000..d2206cf5cf --- /dev/null +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVesselDrawer.cpp @@ -0,0 +1,190 @@ +/*=================================================================== + +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 "mitkPAVesselDrawer.h" +#include "mitkPAVesselProperties.h" + +mitk::pa::VesselDrawer::VesselDrawer() +{ +} + +mitk::pa::VesselDrawer::~VesselDrawer() +{ +} + +void mitk::pa::VesselDrawer::FillInVolumeValues(mitk::pa::InSilicoTissueVolume::Pointer volume, + mitk::pa::VesselProperties::Pointer properties, + double vesselFraction, + int x, int y, int z) +{ + if (vesselFraction <= 1 && vesselFraction >= 0) + { + double backgroundFraction = 1.0 - vesselFraction; + + double absorption = backgroundFraction * + volume->GetAbsorptionVolume()->GetData( + x, y, z) + + vesselFraction * properties->GetAbsorptionCoefficient(); + double scattering = backgroundFraction * + volume->GetScatteringVolume()->GetData( + x, y, z) + + vesselFraction * properties->GetScatteringCoefficient(); + double anisotropy = backgroundFraction * + volume->GetAnisotropyVolume()->GetData( + x, y, z) + + vesselFraction * properties->GetAnisotopyCoefficient(); + + volume->SetVolumeValues(x, y, z, + absorption, + scattering, + anisotropy); + + if (vesselFraction >= 0.5) + { + m_VolumeState->SetData(VolumeState::PARTIAL, x, y, z); + } + } +} + +void mitk::pa::VesselDrawer::InitializeVolumeState(InSilicoTissueVolume::Pointer volume) +{ + Volume::Pointer referenceVolume = volume->GetAbsorptionVolume(); + unsigned int xDim = referenceVolume->GetXDim(); + unsigned int yDim = referenceVolume->GetYDim(); + unsigned int zDim = referenceVolume->GetZDim(); + double* data = new double[xDim*yDim*zDim]; + for (unsigned int idx = 0, length = xDim*yDim*zDim; idx < length; ++idx) + { + data[idx] = VolumeState::UNTOUCHED; + } + + m_VolumeState = Volume::New(data, xDim, yDim, zDim); +} + +void mitk::pa::VesselDrawer::DrawVesselInVolume(Vector::Pointer fromPosition, + VesselProperties::Pointer properties, + InSilicoTissueVolume::Pointer volume) +{ + if (m_VolumeState.IsNull()) + { + InitializeVolumeState(volume); + } + + Vector::Pointer stepSize = properties->GetDirectionVector(); + + 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)) + { + properties->SetRadiusInVoxel(0); + return; + } + + double radius = properties->GetRadiusInVoxel(); + double ceiledRadius = ceil(radius) + 3; + double radiusSquared = radius*radius; + double PI = 3.141; + double totalPoints = 1000; + + 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 - sqrt(xDiff*xDiff + yDiff*yDiff + zDiff*zDiff); + + if (vectorLengthDiff < -2 || + sqrt(xPos*xPos + yPos*yPos + zPos*zPos) > 1) + { + continue; + } + + if (vectorLengthDiff >= 1) + { + volume->SetVolumeValues(x, y, z, + properties->GetAbsorptionCoefficient(), + properties->GetScatteringCoefficient(), + properties->GetAnisotopyCoefficient(), + mitk::pa::InSilicoTissueVolume::SegmentationType::VESSEL); + + m_VolumeState->SetData(VolumeState::FULL, x, y, z); + } + else + { + if (properties->GetDoPartialVolume() && + abs(m_VolumeState->GetData(x, y, z) - VolumeState::UNTOUCHED) < mitk::eps) + { + double 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.0) < -mitk::eps; checkx += 0.1) + for (double checky = 0; (checky - 1.0) < -mitk::eps; checky += 0.1) + for (double checkz = 0; (checkz - 1.0) < -mitk::eps; checkz += 0.1) + { + 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.0; + } + } + + double vesselFraction = (pointsWithin / totalPoints); + + FillInVolumeValues(volume, properties, vesselFraction, x, y, z); + } + } + } +} + +bool mitk::pa::Equal(const VesselDrawer::Pointer leftHandSide, const VesselDrawer::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 (!Equal(leftHandSide->GetVolumeState(), rightHandSide->GetVolumeState(), eps, verbose)) + //{ + // MITK_INFO(verbose) << "Volume state not equal"; + // return false; + //} + + return true; +} diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulatorControls.ui b/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulatorControls.ui index 505c73e918..4187e1e967 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulatorControls.ui +++ b/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulatorControls.ui @@ -1,4725 +1,4728 @@ PASimulatorControls 0 0 437 655 0 0 Ubuntu Qt::NoContextMenu QmitkTemplate :/org.mitk.gui.qt.photoacousticsimulation/resources/icon.xpm:/org.mitk.gui.qt.photoacousticsimulation/resources/icon.xpm 0 0 415 600 Ubuntu 0 Generator 10 10 391 581 0 0 391 581 0 0 0 0 75 true Volume parameters 0 0 0 25 16777215 25 Tissue name: 0 0 0 25 16777215 25 Size x: 0 0 0 25 16777215 25 50 false Spacing: 0 0 0 25 16777215 25 PhotoacousticTissue 0 0 0 25 16777215 25 1 9999 - 140 + 35 0 0 0 25 16777215 25 y: 0 0 0 25 16777215 25 9999 - 200 + 50 0 0 0 25 16777215 25 z: 0 0 0 25 16777215 25 9999 - 200 + 50 0 0 0 25 16777215 25 voxels Qt::Horizontal 40 20 0 0 0 25 16777215 25 - 2 + 5 0.010000000000000 - 0.030000000000000 + 0.120000000000000 0 0 0 25 16777215 25 cm Qt::Horizontal 40 20 0 0 0 25 16777215 25 Randomize: 0 0 0 25 16777215 25 Custom seed: 0 0 0 25 16777215 25 Partial volume effects: Generate batch file output: 0 0 0 25 16777215 25 false true 0 0 0 25 16777215 25 + + true + true 0 0 0 25 16777215 25 true false 0 0 0 25 16777215 25 sigma: 0 0 0 25 16777215 25 0 100.000000000000000 0.010000000000000 2.000000000000000 0 0 0 25 16777215 25 % Qt::Horizontal 40 20 0 0 0 25 16777215 25 999999999 - 170704057 + 1337 Qt::Horizontal 40 20 Qt::Horizontal 40 20 Qt::Horizontal 40 20 Number of volumes to generate: 1 9999999 1 Qt::Horizontal 40 20 Save generated tissue path: 0 0 50 0 50 16777215 open Path to MitkMcxyz binary: 0 0 50 0 50 16777215 open From (cm): 0 0 0 25 16777215 25 -100000.000000000000000 100000.000000000000000 0.100000000000000 -1.800000000000000 To (cm): 0 0 0 25 16777215 25 -10000.000000000000000 10000.000000000000000 0.100000000000000 1.800000000000000 Step: 0 0 0 25 16777215 25 -10000.000000000000000 10000.000000000000000 0.010000000000000 0.120000000000000 Number of Photons (x1000): 0 999999999 1 100000 0 0 Generate Tissue Qt::Vertical 20 40 Tissue 10 10 391 581 0 0 391 581 0 0 0 0 0 0 0 25 16777215 25 11 75 false true false Air Parameters 0 0 0 25 16777215 25 50 false Thickness: 0 0 0 10 16777215 10 11 75 true 0 0 0 25 16777215 25 75 false true false Background Parameters 0 0 0 25 16777215 25 Absorption coefficient: 0 0 0 25 16777215 25 Scattering coefficient: 0 0 0 25 16777215 25 Anisotropy facor: 0 0 0 10 16777215 10 11 75 true 0 0 0 25 16777215 25 75 false true false Skin Parameters 0 0 0 25 16777215 25 Thickness: 0 0 0 25 16777215 25 Absorption coefficient: 0 0 0 25 16777215 25 Scattering coefficient: 0 0 0 25 16777215 25 Anisotropy facor: 0 0 0 10 16777215 10 11 75 true 0 0 0 25 16777215 25 11 75 true 0 0 0 25 16777215 25 999.990000000000009 0.100000000000000 12.000000000000000 0 0 0 25 16777215 25 mm Qt::Horizontal 40 20 0 0 0 10 16777215 10 11 75 true 0 0 0 25 16777215 25 75 true 0 0 0 25 16777215 25 5 0.010000000000000 0.100000000000000 0.100000000000000 0 0 0 25 16777215 25 per cm Qt::Horizontal 40 20 0 0 0 25 16777215 25 5 0.010000000000000 1000.000000000000000 0.500000000000000 15.000000000000000 0 0 0 25 16777215 25 per cm Qt::Horizontal 40 20 0 0 0 25 16777215 25 5 0.010000000000000 1.000000000000000 0.010000000000000 0.900000000000000 Qt::Horizontal 40 20 0 0 0 10 16777215 10 11 75 true 0 0 0 25 16777215 25 75 true 0 0 0 25 16777215 25 0.100000000000000 0.000000000000000 0 0 0 25 16777215 25 mm Qt::Horizontal 40 20 0 0 0 25 16777215 25 5 0.010000000000000 0.100000000000000 3.000000000000000 0 0 0 25 16777215 25 per cm Qt::Horizontal 40 20 0 0 0 25 16777215 25 5 0.010000000000000 1000.000000000000000 0.500000000000000 20.000000000000000 0 0 0 25 16777215 25 per cm Qt::Horizontal 40 20 0 0 0 25 16777215 25 5 0.010000000000000 1.000000000000000 0.010000000000000 0.900000000000000 Qt::Horizontal 40 20 0 0 0 10 16777215 10 11 75 true Qt::Vertical 20 40 Vessels 10 10 391 581 0 0 391 581 0 0 0 0 75 true Vessel Parameters 0 25 16777215 25 The number of bloos vessels to generate Count: Qt::AlignLeading|Qt::AlignLeft|Qt::AlignVCenter 0 25 16777215 25 The radius of the blood vessels in mm Radius: 0 25 16777215 25 the absorption coefficient refers to the number of absorption events per centimeter. Absorption: 0 25 16777215 25 The reduced scattering coefficient. It refers to the amount of scattering events per centimeter. Scattering: 0 25 16777215 25 The anisotropy factor is the probability of a photon to not change its direction after a scattering event. Anisotropy factor: 0 25 16777215 25 The bifurcation frequency determines how often the vessel should bifurcate. The vessel will bifurcate after meandering a mean of the specified amount of voxels. When given a value of 0, the vessel will never bifurcate. Bifurcation frequency: 0 25 16777215 25 The curvedness it a setting to determine how much the vessel is allowed to bend during creation. A value of 0 refers to no bending at all and a value of 5 is the maximum. Curvedness: 0 25 16777215 25 The spawn depth defines the depth range in which the vessels enter the volume. Spawn depth: 0 25 16777215 25 The minimum number of blood vessels 0 1 0 25 16777215 25 to 0 25 16777215 25 The maximum number of blood vessels 1 0 25 16777215 25 Vessels Qt::Horizontal QSizePolicy::Expanding 60 20 0 25 16777215 25 The minimum radius 5 2.250000000000000 0 25 16777215 25 to 0 25 16777215 25 The maximum radius 5 4.050000000000000 0 25 16777215 25 mm Qt::Horizontal QSizePolicy::Expanding 60 20 0 25 16777215 25 The minimum absorption coefficient 5 0.010000000000000 2.000000000000000 0 25 16777215 25 to 0 25 16777215 25 The maximum absorption coefficient 5 0.010000000000000 8.000000000000000 0 25 16777215 25 per cm Qt::Horizontal QSizePolicy::Expanding 60 20 0 25 16777215 25 The minimum scattering 5 0.010000000000000 999.000000000000000 1.000000000000000 15.000000000000000 0 25 16777215 25 to 0 25 16777215 25 The minimum scattering 5 0.010000000000000 999.000000000000000 1.000000000000000 15.000000000000000 0 25 16777215 25 per cm Qt::Horizontal QSizePolicy::Expanding 60 20 0 25 16777215 25 The minimum anisotropy factor 5 0.010000000000000 1.000000000000000 0.100000000000000 0.900000000000000 0 25 16777215 25 to 0 25 16777215 25 The maximum anisotropy factor 5 0.010000000000000 1.000000000000000 0.100000000000000 0.900000000000000 Qt::Horizontal QSizePolicy::Expanding 60 20 0 25 16777215 25 The bifurcation frequency determines how often the vessel should bifurcate. The vessel will bifurcate after meandering a mean of the specified amount of voxels. When given a value of 0, the vessel will never bifurcate. 0 1.000000000000000 999999999.000000000000000 5.000000000000000 50.000000000000000 0 25 16777215 25 voxels Qt::Horizontal QSizePolicy::Expanding 60 20 0 25 16777215 25 The minimal curvedness of the vessel. A value of 0 refers to no bending at all and a value of 5 is the maximum. 5.000000000000000 0.000000000000000 0 25 16777215 25 to 0 25 16777215 25 The maximal curvedness of the vessel. A value of 0 refers to no bending at all and a value of 5 is the maximum. 5.000000000000000 0.010000000000000 0.200000000000000 0 25 16777215 25 A.U. Qt::Horizontal QSizePolicy::Expanding 60 20 0 25 16777215 25 The minimum spawn depth 5 0.010000000000000 2.200000000000000 0 25 16777215 25 to 0 25 16777215 25 the maximum spawn depth 5 0.010000000000000 4.200000000000000 0 25 16777215 25 cm Qt::Horizontal QSizePolicy::Expanding 60 20 Qt::Vertical 20 40 Monte Carlo 10 10 391 581 0 0 391 581 0 0 0 0 75 true Monte Carlo Parameters 0 25 16777215 25 50 false General: 0 0 0 25 16777215 25 Mcflag: 0 0 0 25 16777215 25 Launchflag: 0 0 0 25 16777215 25 Boundaryflag: 0 0 0 25 16777215 25 1 4 Qt::Horizontal 40 20 0 0 0 25 16777215 25 0 0 Qt::Horizontal 40 20 0 0 0 25 16777215 25 1 2 Qt::Horizontal 40 20 0 25 16777215 25 50 false Initial launch point: 0 25 16777215 25 x 0 25 16777215 25 4 1000000.000000000000000 0.000000000000000 0 25 16777215 25 y 0 25 16777215 25 4 1000000.000000000000000 0.000000000000000 0 25 16777215 25 z 0 25 16777215 25 4 1000000.000000000000000 0.000000000000000 Qt::Horizontal 40 20 0 25 16777215 25 50 false Focus point: 0 25 16777215 25 x 0 25 16777215 25 4 1000000.000000000000000 0.000000000000000 0 25 16777215 25 y 0 25 16777215 25 4 1000000.000000000000000 0.000000000000000 0 25 16777215 25 z 0 25 16777215 25 4 1000000.000000000000000 0.000000000000000 Qt::Horizontal 40 20 0 25 16777215 25 50 false Trajectory vector: 0 25 16777215 25 x 0 25 16777215 25 4 1000000.000000000000000 0.000000000000000 0 25 16777215 25 y 0 25 16777215 25 4 1000000.000000000000000 0.342000000000000 0 25 16777215 25 z 0 25 16777215 25 4 1000000.000000000000000 0.939700000000000 Qt::Horizontal 40 20 radius: waist: 4 1000.000000000000000 0.500000000000000 Qt::Horizontal 40 20 4 1000.000000000000000 0.010000000000000 Qt::Horizontal 40 20 Qt::Vertical 20 40 Wavelength 10 10 391 581 0 0 391 581 Qt::NoContextMenu 0 0 0 0 75 true Adjust physical properties by wavelength false 16777215 250 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0//EN" "http://www.w3.org/TR/REC-html40/strict.dtd"> <html><head><meta name="qrichtext" content="1" /><style type="text/css"> p, li { white-space: pre-wrap; } </style></head><body style=" font-family:'Ubuntu'; font-size:7.8pt; font-weight:400; font-style:normal;"> <p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-size:11pt;">This widget enables the adjustment of the physical properties of the tissue according to a selected wavelength of the light and the oxygen saturation of the blood.</span></p> <p style="-qt-paragraph-type:empty; margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px; font-size:11pt;"><br /></p> <p style=" margin-top:0px; margin-bottom:0px; margin-left:0px; margin-right:0px; -qt-block-indent:0; text-indent:0px;"><span style=" font-size:11pt;">The algorithm and measurements were provided by the publication </span><span style=" font-size:11pt; font-weight:600;">Optical properties of biological tissues: a review </span><span style=" font-size:11pt;">by Steve L. Jacques.</span></p></body></html> 0 0 0 25 16777215 25 Wavelength: 0 0 0 25 16777215 25 Vessel oxygen saturation: 0 0 0 25 16777215 25 0 300.000000000000000 1000.000000000000000 650.000000000000000 Qt::Horizontal 40 20 0 0 0 25 16777215 25 0 100.000000000000000 75.000000000000000 0 0 0 25 16777215 25 % Qt::Horizontal 40 20 Adjust tissue properties Qt::Vertical 20 40