diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp index 705f262b2b..f9ccad987a 100644 --- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp +++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp @@ -1,162 +1,175 @@ /*=================================================================== 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 "mitkPhotoacousticTissueGenerator.h" #include "mitkPhotoacousticSmartVector.h" #include "mitkPhotoacousticFatTissueGenerator.h" mitk::PhotoacousticVolume::Pointer mitk::PhotoacousticTissueGenerator::GenerateInSilicoData(int xDim, int yDim, int zDim, double spacing, double basicAbsorption, double basicScattering, double basicAnisotropy, double airThickness, double skinAbsorption, double skinScattering, double skinAnisotropy, double skinThickness, double fatAbsorption, double fatScattering, double fatAnisotropy, double fatPercentage, bool randomizePhysicalProperties, double randomPercentage, int minNumberOfVessels, int maxNumberOfVessels, double minBending, double maxBending, double minAbs, double maxAbs, double minRadius, double maxRadius, double mcflag, double launchflag, double boundaryflag, double launchPointX, double launchPointY, double launchPointZ, double focusPointX, double focusPointY, double focusPointZ, double trajectoryVectorX, double trajectoryVectorY, double trajectoryVectorZ, double radius, double waist, double sigma, bool doGauss, int bifurcationFrequency, - double vesselScattering, double vesselAnisotropy, + double vesselScatteringMinimum, double vesselScatteringMaximum, + double vesselAnisotropyMinimum, double vesselAnisotropyMaximum, bool useRngSeed, long rngSeed, mitk::PhotoacousticVessel::CalculateNewVesselPositionCallback calculateNewVesselPositionCallback) { MITK_DEBUG << "Initializing Empty Volume"; + mitk::PhotoacousticVolume::Pointer generatedVolume = mitk::PhotoacousticVolume::New(); + generatedVolume->Initialize(xDim, yDim, zDim, spacing, basicAbsorption, basicScattering, basicAnisotropy, mcflag, launchflag, boundaryflag, launchPointX, launchPointY, launchPointZ, focusPointX, focusPointY, focusPointZ, trajectoryVectorX, trajectoryVectorY, trajectoryVectorZ, radius, waist, sigma, doGauss, airThickness, skinThickness, skinAbsorption, skinScattering, skinAnisotropy); + const double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2; + std::mt19937 randomNumberGenerator; std::random_device randomDevice; if(useRngSeed) { randomNumberGenerator.seed(rngSeed); } else { if(randomDevice.entropy()>0.1) { randomNumberGenerator.seed(randomDevice()); } else { randomNumberGenerator.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count()); } } mitk::PhotoacousticFatTissueGenerator::CreateFatGlobules(generatedVolume, fatPercentage, fatAbsorption, fatScattering, fatAnisotropy, &randomNumberGenerator); std::uniform_int_distribution randomNumVesselDistribution(minNumberOfVessels, maxNumberOfVessels); std::uniform_real_distribution randomBendingDistribution(minBending, maxBending); std::uniform_real_distribution randomAbsorptionDistribution(minAbs, maxAbs); std::uniform_real_distribution randomRadiusDistribution(minRadius, maxRadius); - std::uniform_real_distribution randomSign(-1, 1); + std::uniform_real_distribution randomScatteringDistribution(vesselScatteringMinimum, vesselScatteringMaximum); + std::uniform_real_distribution randomAnisotropyDistribution(vesselAnisotropyMinimum, vesselAnisotropyMaximum); + std::uniform_int_distribution borderTypeDistribution(0, 3); int numberOfBloodVessels = randomNumVesselDistribution(randomNumberGenerator); generatedVolume->AddIntProperty("numberOfVessels", numberOfBloodVessels); generatedVolume->AddIntProperty("bifurcationFrequency", bifurcationFrequency); + int airvoxels = airThickness / spacing / 10; + MITK_DEBUG << "Simulating vessel structures"; for(int vesselNumber = 0; vesselNumber < numberOfBloodVessels; vesselNumber++) { mitk::PhotoacousticSmartVector::Pointer initialPosition = mitk::PhotoacousticSmartVector::New(); mitk::PhotoacousticSmartVector::Pointer initialDirection = mitk::PhotoacousticSmartVector::New(); double initialRadius = randomRadiusDistribution(randomNumberGenerator) / spacing / 10; std::stringstream radiusString; radiusString <<"vessel_"<AddDoubleProperty(radiusString.str(), initialRadius); double absorptionCoefficient = randomAbsorptionDistribution(randomNumberGenerator); std::stringstream absorptionString; absorptionString <<"vessel_"<AddDoubleProperty(absorptionString.str(), absorptionCoefficient); double bendingFactor = randomBendingDistribution(randomNumberGenerator); std::stringstream bendingString; bendingString <<"vessel_"<AddDoubleProperty(bendingString.str(), bendingFactor); - bool stickXBorder = randomSign(randomNumberGenerator) < 0; - bool stickLeft = randomSign(randomNumberGenerator) < 0; - - if(stickXBorder) - { - if(stickLeft) - { - initialPosition->Randomize(0, 0, initialRadius, yDim-initialRadius, 3+initialRadius, zDim-initialRadius, &randomNumberGenerator); - initialDirection->Randomize(1, 2, -DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM, -DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM,&randomNumberGenerator); - } - else - { - initialPosition->Randomize(xDim, xDim, initialRadius, yDim - initialRadius, 3+initialRadius, zDim, &randomNumberGenerator); - initialDirection->Randomize(-2, -1, -DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM, -DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM, &randomNumberGenerator); - } - } - else - { - if(stickLeft) - { - initialPosition->Randomize(initialRadius, xDim - initialRadius, 0, 0, 3+initialRadius, zDim - initialRadius, &randomNumberGenerator); - initialDirection->Randomize(-DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM, 1, 2, -DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM, &randomNumberGenerator); - } - else - { - initialPosition->Randomize(initialRadius, xDim - initialRadius, yDim, yDim, 3+initialRadius, zDim - initialRadius, &randomNumberGenerator); - initialDirection->Randomize(-DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM, -2, -1, -DIRECTION_VECTOR_WIGGLE_ROOM, DIRECTION_VECTOR_WIGGLE_ROOM, &randomNumberGenerator); - } - } - - mitk::PhotoacousticVesselTree::Pointer vesselStructure = mitk::PhotoacousticVesselTree::New(); - + double vesselScattering = randomScatteringDistribution(randomNumberGenerator); std::stringstream scatteringString; scatteringString <<"vessel_"<AddDoubleProperty(scatteringString.str(), vesselScattering); + double vesselAnisotropy= randomAnisotropyDistribution(randomNumberGenerator); std::stringstream anisotropyString; anisotropyString <<"vessel_"<AddDoubleProperty(anisotropyString.str(), vesselAnisotropy); + /*The vessel tree shall start at one of the 4 sides of the volume. + * The vessels will always be completely contained in the volume + * when starting to meander. + * They will meander in a direction perpendicular to the + * plane they started from (within the limits of the + * DIRECTION_VECTOR_INITIAL_VARIANCE) + */ + int borderType = borderTypeDistribution(randomNumberGenerator); + switch(borderType) + { + case 0: + initialPosition->Randomize(0, 0, initialRadius, yDim-initialRadius, airvoxels+initialRadius, zDim-initialRadius, &randomNumberGenerator); + initialDirection->Randomize(1, 2, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, + -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE,&randomNumberGenerator); + break; + case 1: + initialPosition->Randomize(xDim, xDim, initialRadius, yDim - initialRadius, airvoxels+initialRadius, zDim, &randomNumberGenerator); + initialDirection->Randomize(-2, -1, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, + -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); + break; + case 2: + initialPosition->Randomize(initialRadius, xDim - initialRadius, 0, 0, airvoxels+initialRadius, zDim - initialRadius, &randomNumberGenerator); + initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, 1, 2, + -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); + break; + case 3: + initialPosition->Randomize(initialRadius, xDim - initialRadius, yDim, yDim, airvoxels+initialRadius, zDim - initialRadius, &randomNumberGenerator); + initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -2, -1, + -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); + break; + } + + mitk::PhotoacousticVesselTree::Pointer vesselTree = mitk::PhotoacousticVesselTree::New(); + MITK_DEBUG << "Initializing vessel structure"; - vesselStructure->Initialize(initialPosition, initialDirection, initialRadius, absorptionCoefficient, vesselScattering, vesselAnisotropy, bifurcationFrequency, &randomNumberGenerator); + vesselTree->Initialize(initialPosition, initialDirection, initialRadius, absorptionCoefficient, vesselScattering, vesselAnisotropy, bifurcationFrequency, &randomNumberGenerator); MITK_DEBUG << "Simulating vessel structure"; - while(!vesselStructure->IsFinished()) + while(!vesselTree->IsFinished()) { - vesselStructure->Step(generatedVolume, calculateNewVesselPositionCallback, bendingFactor, &randomNumberGenerator); + vesselTree->Step(generatedVolume, calculateNewVesselPositionCallback, bendingFactor, &randomNumberGenerator); } } generatedVolume->GaussianBlur3D(); generatedVolume->FinalizeVolume(randomizePhysicalProperties, randomPercentage, useRngSeed, rngSeed); return generatedVolume; } diff --git a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.h index 5b9f036d73..c6221df255 100644 --- a/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.h +++ b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.h @@ -1,114 +1,113 @@ /*=================================================================== 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 mitkPhotoacousticTissueGenerator_h #define mitkPhotoacousticTissueGenerator_h #include #include #include #include #include #include "mitkPhotoacousticVesselTree.h" #include "mitkPhotoacousticVolume.h" #include "mitkCommon.h" namespace mitk { class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticTissueGenerator final { public: /** * @brief GenerateInSilicoData This method will return a PhotoacousticVolume created in terms of the given parameters. * @param xDim x dimension of the generated volume * @param yDim y dimension of the generated volume * @param zDim z dimension of the generated volume * @param spacing voxel spacing. Every dimension is equally spaced. * @param basicAbsorption absorption coefficient for standard tissue * @param basicScattering scattering coefficient for standard tissue * @param basicAnisotropy anisotropy factor for standard tissue * @param airThickness * @param skinAbsorption * @param skinScattering * @param skinAnisotropy * @param skinThickness * @param fatAbsorption * @param fatScattering * @param fatAnisotropy * @param fatPercentage, * @param minNumberOfVessels minimum number of vessels generated in volume * @param maxNumberOfVessels maximum number of vessels generated in volume * @param minBending minimum bending coefficient (0 = not bending; >1 = bending a lot) * @param maxBending maximum bending coefficient * @param minAbs minimum absorption coefficient for blood vessel * @param maxAbs maximum absorption coefficient for blood vessel * @param minRadius minimum vessel radius * @param maxRadius maximum vessel radius * @param mcflag the montecarlo flag (chosen probe design) * @param launchflag * @param boundaryflag * @param launchPointX * @param launchPointY * @param launchPointZ * @param focusPointX * @param focusPointY * @param focusPointZ * @param trajectoryVectorX * @param trajectoryVectorY * @param trajectoryVectorZ * @param radius * @param waist * @param sigma * @param doGauss * @param bifurcationFrequency * @param calculateNewVesselPositionCallback * @return */ static mitk::PhotoacousticVolume::Pointer GenerateInSilicoData(int xDim, int yDim, int zDim, double spacing, double basicAbsorption, double basicScattering, double basicAnisotropy, double airThickness, double skinAbsorption, double skinScattering, double skinAnisotropy, double skinThickness, double fatAbsorption, double fatScattering, double fatAnisotropy, double fatPercentage, bool randomizePhysicalproperties, double randomPercentage, int minNumberOfVessels, int maxNumberOfVessels, double minBending, double maxBending, double minAbs, double maxAbs, double minRadius, double maxRadius, double mcflag, double launchflag, double boundaryflag, double launchPointX, double launchPointY, double launchPointZ, double focusPointX, double focusPointY, double focusPointZ, double trajectoryVectorX, double trajectoryVectorY, double trajectoryVectorZ, double radius, double waist, double sigma, bool doGauss, int bifurcationFrequency, - double vesselScattering, double vesselAnisotropy, + double vesselScatteringMinimum, double vesselScatteringMaximum, + double vesselAnisotropyMinimum, double vesselAnisotropyMaximum, bool useRngSeed, long rngSeed, mitk::PhotoacousticVessel::CalculateNewVesselPositionCallback calculateNewVesselPositionCallback); private: PhotoacousticTissueGenerator(); virtual ~PhotoacousticTissueGenerator(); - const double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2; - }; } #endif