diff --git a/Modules/PhotoacousticsLib/CMakeLists.txt b/Modules/PhotoacousticsLib/CMakeLists.txt index 0ae41f8f5c..c11c71a190 100644 --- a/Modules/PhotoacousticsLib/CMakeLists.txt +++ b/Modules/PhotoacousticsLib/CMakeLists.txt @@ -1,15 +1,16 @@ MITK_CREATE_MODULE( INCLUDE_DIRS PUBLIC include INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} DEPENDS PUBLIC MitkAlgorithmsExt tinyxml PACKAGE_DEPENDS tinyxml PUBLIC ITK ) add_subdirectory(MitkMCxyz) add_subdirectory(MitkTissueBatchGenerator) add_subdirectory(MitkPAPhantomGenerator) +add_subdirectory(MitkOxygenationTissueBatchGenerator) add_subdirectory(test) diff --git a/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/CMakeLists.txt b/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/CMakeLists.txt new file mode 100644 index 0000000000..7b3c2e6b39 --- /dev/null +++ b/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/CMakeLists.txt @@ -0,0 +1,11 @@ +OPTION(BUILD_PhotoacousticTissueBatchGenerator "Build MiniApp for batch generating of photoacoustic tissue in different wavelengths" ON) + +IF(BUILD_PhotoacousticTissueBatchGenerator) + PROJECT( MitkOxygenationTissueBatchGenerator ) + mitk_create_executable(OxygenationBatchGenerator + DEPENDS MitkCommandLine MitkCore MitkPhotoacousticsLib + PACKAGE_DEPENDS + CPP_FILES OxygenationBatchGenerator.cpp) + + install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) + ENDIF() diff --git a/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/OxygenationBatchGenerator.cpp b/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/OxygenationBatchGenerator.cpp new file mode 100644 index 0000000000..3cb80bdf27 --- /dev/null +++ b/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/OxygenationBatchGenerator.cpp @@ -0,0 +1,258 @@ +/*=================================================================== + +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 +#include +#include +#include +#include +#include +#include + +#include + +using namespace mitk::pa; + +TissueGeneratorParameters::Pointer CreateOxygenationParameters(int wavelength, double oxygenation, long rng) +{ + auto propertyCalculator = PropertyCalculator::New(); + auto background_properties = propertyCalculator->CalculatePropertyForSpecificWavelength( + PropertyCalculator::TissueType::STANDARD_TISSUE, wavelength, oxygenation); + auto returnParameters = TissueGeneratorParameters::New(); + returnParameters->SetAirThicknessInMillimeters(12); + returnParameters->SetMinBackgroundAbsorption(background_properties.mua); + returnParameters->SetMaxBackgroundAbsorption(background_properties.mua); + returnParameters->SetBackgroundAnisotropy(background_properties.g); + returnParameters->SetBackgroundScattering(background_properties.mus); + returnParameters->SetCalculateNewVesselPositionCallback(&VesselMeanderStrategy::CalculateRandomlyDivergingPosition); + returnParameters->SetDoPartialVolume(false); + returnParameters->SetMinNumberOfVessels(2); + returnParameters->SetMaxNumberOfVessels(10); + //---These are calculated for each vessel individually + returnParameters->SetMinVesselAbsorption(0); + returnParameters->SetMaxVesselAbsorption(0); + returnParameters->SetMinVesselAnisotropy(0); + returnParameters->SetMaxVesselAnisotropy(0); + returnParameters->SetMinVesselScattering(0); + returnParameters->SetMaxVesselScattering(0); + //--- + returnParameters->SetMinVesselBending(0.1); + returnParameters->SetMaxVesselBending(0.3); + returnParameters->SetMinVesselRadiusInMillimeters(0.5); + returnParameters->SetMaxVesselRadiusInMillimeters(4); + returnParameters->SetMinVesselZOrigin(2.2); + returnParameters->SetMaxVesselZOrigin(4); + returnParameters->SetVesselBifurcationFrequency(50); + returnParameters->SetRandomizePhysicalProperties(false); + returnParameters->SetSkinThicknessInMillimeters(0); + returnParameters->SetUseRngSeed(true); + returnParameters->SetRngSeed(rng); + returnParameters->SetVoxelSpacingInCentimeters(0.06); + returnParameters->SetXDim(70); + returnParameters->SetYDim(50); + returnParameters->SetZDim(90); + returnParameters->SetRandomizeVesselOpticalParams(true); + returnParameters->SetTargetWavelength(wavelength); + return returnParameters; +} + +struct InputParameters +{ + std::string saveFolderPath; + std::string identifyer; + std::string exePath; + std::string probePath; + bool verbose; +}; + +InputParameters parseInput(int argc, char* argv[]) +{ + MITK_INFO << "Paring arguments..."; + mitkCommandLineParser parser; + // set general information + parser.setCategory("MITK-Photoacoustics"); + parser.setTitle("Mitk Oxygenation Tissue Batch Generator"); + parser.setDescription("Creates in silico tissue in batch processing and automatically calculates fluence values for the central slice of the volume."); + parser.setContributor("Computer Assisted Medical Interventions, DKFZ"); + + // how should arguments be prefixed + parser.setArgumentPrefix("--", "-"); + // add each argument, unless specified otherwise each argument is optional + // see mitkCommandLineParser::addArgument for more information + parser.beginGroup("Required parameters"); + parser.addArgument( + "savePath", "s", mitkCommandLineParser::InputDirectory, + "Input save folder (directory)", "input save folder", + us::Any(), false); + parser.addArgument( + "mitkMcxyz", "m", mitkCommandLineParser::OutputFile, + "MitkMcxyz binary (file)", "path to the MitkMcxyz binary", + us::Any(), false); + parser.endGroup(); + parser.beginGroup("Optional parameters"); + parser.addArgument( + "probe", "p", mitkCommandLineParser::OutputFile, + "xml probe file (file)", "file to the definition of the used probe (*.xml)", + us::Any()); + parser.addArgument( + "verbose", "v", mitkCommandLineParser::Bool, + "Verbose Output", "Whether to produce verbose, or rather debug output"); + parser.addArgument( + "identifyer", "i", mitkCommandLineParser::String, + "Generator identifyer (string)", "A unique identifyer for the calculation instance"); + + InputParameters input; + + std::map parsedArgs = parser.parseArguments(argc, argv); + if (parsedArgs.size() == 0) + exit(-1); + + if (parsedArgs.count("verbose")) + { + MITK_INFO << "verbose"; + input.verbose = us::any_cast(parsedArgs["verbose"]); + } + else + { + input.verbose = false; + } + + if (parsedArgs.count("savePath")) + { + MITK_INFO << "savePath"; + input.saveFolderPath = us::any_cast(parsedArgs["savePath"]); + } + + if (parsedArgs.count("mitkMcxyz")) + { + MITK_INFO << "mitkMcxyz"; + input.exePath = us::any_cast(parsedArgs["mitkMcxyz"]); + } + + if (parsedArgs.count("probe")) + { + MITK_INFO << "probe"; + input.probePath = us::any_cast(parsedArgs["probe"]); + } + + if (parsedArgs.count("identifyer")) + { + MITK_INFO << "identifyer"; + input.identifyer = us::any_cast(parsedArgs["identifyer"]); + } + else + { + MITK_INFO << "generating identifyer"; + auto uid = mitk::UIDGenerator("", 8); + input.identifyer = uid.GetUID(); + } + MITK_INFO << "Paring arguments...[Done]"; + return input; +} + +void simulate(int wavelength, double oxygenation, long rng_seed, int iterationNumber, InputParameters input) +{ + auto parameters = CreateOxygenationParameters(wavelength, oxygenation, rng_seed); + MITK_INFO(input.verbose) << "Generating tissue.."; + auto resultTissue = InSilicoTissueGenerator::GenerateInSilicoData(parameters); + MITK_INFO(input.verbose) << "Generating tissue..[Done]"; + + auto inputfolder = std::string(input.saveFolderPath + "input/"); + auto outputfolder = std::string(input.saveFolderPath + "output/"); + if (!itksys::SystemTools::FileIsDirectory(inputfolder)) + { + itksys::SystemTools::MakeDirectory(inputfolder); + } + if (!itksys::SystemTools::FileIsDirectory(outputfolder)) + { + itksys::SystemTools::MakeDirectory(outputfolder); + } + + std::string savePath = input.saveFolderPath + "input/BaselineHB_" + input.identifyer + + "_" + std::to_string(iterationNumber) + "_oxy" + std::to_string(oxygenation); + std::string outputPath = input.saveFolderPath + "output/BaselineHB_" + input.identifyer + + "_" + std::to_string(iterationNumber) + "_oxy" + std::to_string(oxygenation); + + if (!itksys::SystemTools::FileIsDirectory(outputPath)) + { + itksys::SystemTools::MakeDirectory(outputPath); + } + if (!itksys::SystemTools::FileIsDirectory(savePath)) + { + itksys::SystemTools::MakeDirectory(savePath); + } + + savePath = savePath + "/wl" + std::to_string(wavelength) + ".nrrd"; + mitk::IOUtil::Save(resultTissue->ConvertToMitkImage(), savePath); + + outputPath = outputPath + "/wl" + std::to_string(wavelength); + + MITK_INFO(input.verbose) << "Simulating fluence.."; + double yo=0; + std::string yo_string = std::to_string(round(yo*100)/100.0); + int result = -4; + if(!input.probePath.empty()) + result = std::system(std::string(input.exePath + " -i " + savePath + " -o " + + (outputPath + ".nrrd") + + " -yo " + yo_string + " -p " + input.probePath + + " -n 10000000").c_str()); + else + result = std::system(std::string(input.exePath + " -i " + savePath + " -o " + + (outputPath + "_oxy" + std::to_string((int) (oxygenation*100)) + + "_wl" + std::to_string(wavelength) + ".nrrd") + + " -yo " + yo_string + " -n 10000000").c_str()); + MITK_INFO << "wl: " << wavelength << ": " << result; + + MITK_INFO(input.verbose) << "Simulating fluence..[Done]"; +} + +int main(int argc, char * argv[]) +{ + auto input = parseInput(argc, argv); + unsigned int iterationNumber = 0; + + std::mt19937 randomNumberGenerator; + std::random_device randomDevice; + + + 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()); + } + + std::uniform_real_distribution randomOxygenation(0, 1); + std::uniform_int_distribution randomLong(0, std::numeric_limits::max()); + + while (true) + { + long rng_seed = randomLong(randomNumberGenerator); + double oxygenation = randomOxygenation(randomNumberGenerator); + + for(int wavelength = 700; wavelength < 960; wavelength+=10) + { + simulate(wavelength, oxygenation, rng_seed, iterationNumber, input); + } + + iterationNumber++; + } +} diff --git a/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/files.cmake b/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/files.cmake new file mode 100644 index 0000000000..9282c291e3 --- /dev/null +++ b/Modules/PhotoacousticsLib/MitkOxygenationTissueBatchGenerator/files.cmake @@ -0,0 +1,3 @@ +set(CPP_FILES + OxygenationBatchGenerator.cpp +) diff --git a/Modules/PhotoacousticsLib/include/mitkPATissueGenerator.h b/Modules/PhotoacousticsLib/include/mitkPATissueGenerator.h index 92f7ab0610..bedce662b4 100644 --- a/Modules/PhotoacousticsLib/include/mitkPATissueGenerator.h +++ b/Modules/PhotoacousticsLib/include/mitkPATissueGenerator.h @@ -1,52 +1,53 @@ /*=================================================================== 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 #include "mitkPAVesselTree.h" #include "mitkPAInSilicoTissueVolume.h" #include "mitkCommon.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT InSilicoTissueGenerator final { public: /** * @brief GenerateInSilicoData This method will return a InSilicoTissueVolume created in terms of the given parameters. * @param parameters * @return */ static InSilicoTissueVolume::Pointer GenerateInSilicoData(TissueGeneratorParameters::Pointer parameters); private: InSilicoTissueGenerator(); virtual ~InSilicoTissueGenerator(); + }; } } #endif diff --git a/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h b/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h index 2ab5fcf444..2cddd5380a 100644 --- a/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h +++ b/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h @@ -1,217 +1,226 @@ /*=================================================================== 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 MITKPHOTOACOUSTICTISSUEGENERATORPARAMETERS_H #define MITKPHOTOACOUSTICTISSUEGENERATORPARAMETERS_H #include #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT TissueGeneratorParameters : public itk::Object { public: mitkClassMacroItkParent(TissueGeneratorParameters, itk::Object) itkFactorylessNewMacro(Self) /** * Callback function definition of a VesselMeanderStrategy */ typedef void (VesselMeanderStrategy::*CalculateNewVesselPositionCallback) (Vector::Pointer, Vector::Pointer, double, std::mt19937*); itkGetMacro(XDim, int) itkGetMacro(YDim, int) itkGetMacro(ZDim, int) itkGetMacro(VoxelSpacingInCentimeters, double) itkGetMacro(DoPartialVolume, bool) itkGetMacro(UseRngSeed, bool) itkGetMacro(RngSeed, long) itkGetMacro(RandomizePhysicalProperties, bool) itkGetMacro(RandomizePhysicalPropertiesPercentage, double) itkGetMacro(ForceVesselsMoveAlongYDirection, bool) itkGetMacro(MinBackgroundAbsorption, double) itkGetMacro(MaxBackgroundAbsorption, double) itkGetMacro(BackgroundScattering, double) itkGetMacro(BackgroundAnisotropy, double) itkGetMacro(AirAbsorption, double) itkGetMacro(AirScattering, double) itkGetMacro(AirAnisotropy, double) itkGetMacro(AirThicknessInMillimeters, double) itkGetMacro(SkinAbsorption, double) itkGetMacro(SkinScattering, double) itkGetMacro(SkinAnisotropy, double) itkGetMacro(SkinThicknessInMillimeters, double) itkGetMacro(CalculateNewVesselPositionCallback, CalculateNewVesselPositionCallback) itkGetMacro(MinNumberOfVessels, int) itkGetMacro(MaxNumberOfVessels, int) itkGetMacro(MinVesselBending, double) itkGetMacro(MaxVesselBending, double) itkGetMacro(MinVesselAbsorption, double) itkGetMacro(MaxVesselAbsorption, double) itkGetMacro(MinVesselRadiusInMillimeters, double) itkGetMacro(MaxVesselRadiusInMillimeters, double) itkGetMacro(VesselBifurcationFrequency, int) itkGetMacro(MinVesselScattering, double) itkGetMacro(MaxVesselScattering, double) itkGetMacro(MinVesselAnisotropy, double) itkGetMacro(MaxVesselAnisotropy, double) itkGetMacro(MinVesselZOrigin, double) itkGetMacro(MaxVesselZOrigin, double) itkGetMacro(MCflag, double) itkGetMacro(MCLaunchflag, double) itkGetMacro(MCBoundaryflag, double) itkGetMacro(MCLaunchPointX, double) itkGetMacro(MCLaunchPointY, double) itkGetMacro(MCLaunchPointZ, double) itkGetMacro(MCFocusPointX, double) itkGetMacro(MCFocusPointY, double) itkGetMacro(MCFocusPointZ, double) itkGetMacro(MCTrajectoryVectorX, double) itkGetMacro(MCTrajectoryVectorY, double) itkGetMacro(MCTrajectoryVectorZ, double) itkGetMacro(MCRadius, double) itkGetMacro(MCWaist, double) + itkGetMacro(RandomizeVesselOpticalParams, bool) + itkGetMacro(TargetWavelength, int) + itkSetMacro(XDim, int) itkSetMacro(YDim, int) itkSetMacro(ZDim, int) itkSetMacro(VoxelSpacingInCentimeters, double) itkSetMacro(DoPartialVolume, bool) itkSetMacro(UseRngSeed, bool) itkSetMacro(RngSeed, long) itkSetMacro(RandomizePhysicalProperties, bool) itkSetMacro(RandomizePhysicalPropertiesPercentage, double) itkSetMacro(ForceVesselsMoveAlongYDirection, bool) itkSetMacro(MinBackgroundAbsorption, double) itkSetMacro(MaxBackgroundAbsorption, double) itkSetMacro(BackgroundScattering, double) itkSetMacro(BackgroundAnisotropy, double) itkSetMacro(AirAbsorption, double) itkSetMacro(AirScattering, double) itkSetMacro(AirAnisotropy, double) itkSetMacro(AirThicknessInMillimeters, double) itkSetMacro(SkinAbsorption, double) itkSetMacro(SkinScattering, double) itkSetMacro(SkinAnisotropy, double) itkSetMacro(SkinThicknessInMillimeters, double) itkSetMacro(CalculateNewVesselPositionCallback, CalculateNewVesselPositionCallback) itkSetMacro(MinNumberOfVessels, int) itkSetMacro(MaxNumberOfVessels, int) itkSetMacro(MinVesselBending, double) itkSetMacro(MaxVesselBending, double) itkSetMacro(MinVesselAbsorption, double) itkSetMacro(MaxVesselAbsorption, double) itkSetMacro(MinVesselRadiusInMillimeters, double) itkSetMacro(MaxVesselRadiusInMillimeters, double) itkSetMacro(VesselBifurcationFrequency, int) itkSetMacro(MinVesselScattering, double) itkSetMacro(MaxVesselScattering, double) itkSetMacro(MinVesselAnisotropy, double) itkSetMacro(MaxVesselAnisotropy, double) itkSetMacro(MinVesselZOrigin, double) itkSetMacro(MaxVesselZOrigin, double) itkSetMacro(MCflag, double) itkSetMacro(MCLaunchflag, double) itkSetMacro(MCBoundaryflag, double) itkSetMacro(MCLaunchPointX, double) itkSetMacro(MCLaunchPointY, double) itkSetMacro(MCLaunchPointZ, double) itkSetMacro(MCFocusPointX, double) itkSetMacro(MCFocusPointY, double) itkSetMacro(MCFocusPointZ, double) itkSetMacro(MCTrajectoryVectorX, double) itkSetMacro(MCTrajectoryVectorY, double) itkSetMacro(MCTrajectoryVectorZ, double) itkSetMacro(MCRadius, double) itkSetMacro(MCWaist, double) + itkSetMacro(RandomizeVesselOpticalParams, bool) + itkSetMacro(TargetWavelength, int) + protected: TissueGeneratorParameters(); ~TissueGeneratorParameters(); private: int m_XDim; int m_YDim; int m_ZDim; double m_VoxelSpacingInCentimeters; bool m_DoPartialVolume; bool m_UseRngSeed; long m_RngSeed; bool m_RandomizePhysicalProperties; double m_RandomizePhysicalPropertiesPercentage; bool m_ForceVesselsMoveAlongYDirection; double m_MinBackgroundAbsorption; double m_MaxBackgroundAbsorption; double m_BackgroundScattering; double m_BackgroundAnisotropy; double m_AirAbsorption; double m_AirScattering; double m_AirAnisotropy; double m_AirThicknessInMillimeters; double m_SkinAbsorption; double m_SkinScattering; double m_SkinAnisotropy; double m_SkinThicknessInMillimeters; + bool m_RandomizeVesselOpticalParams; + int m_TargetWavelength; + CalculateNewVesselPositionCallback m_CalculateNewVesselPositionCallback; int m_MinNumberOfVessels; int m_MaxNumberOfVessels; double m_MinVesselBending; double m_MaxVesselBending; double m_MinVesselAbsorption; double m_MaxVesselAbsorption; double m_MinVesselRadiusInMillimeters; double m_MaxVesselRadiusInMillimeters; int m_VesselBifurcationFrequency; double m_MinVesselScattering; double m_MaxVesselScattering; double m_MinVesselAnisotropy; double m_MaxVesselAnisotropy; double m_MinVesselZOrigin; double m_MaxVesselZOrigin; double m_MCflag; double m_MCLaunchflag; double m_MCBoundaryflag; double m_MCLaunchPointX; double m_MCLaunchPointY; double m_MCLaunchPointZ; double m_MCFocusPointX; double m_MCFocusPointY; double m_MCFocusPointZ; double m_MCTrajectoryVectorX; double m_MCTrajectoryVectorY; double m_MCTrajectoryVectorZ; double m_MCRadius; double m_MCWaist; }; } } #endif // MITKPHOTOACOUSTICTISSUEGENERATORPARAMETERS_H diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp index cbe23af73d..4785adc8be 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp @@ -1,199 +1,230 @@ /*=================================================================== 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 "mitkPATissueGenerator.h" #include "mitkPAVector.h" #include "mitkPAVolumeManipulator.h" +#include "mitkPAPropertyCalculator.h" mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::InSilicoTissueGenerator::GenerateInSilicoData( TissueGeneratorParameters::Pointer parameters) { MITK_DEBUG << "Initializing Empty Volume"; const double RESAMPLING_FACTOR = 2; if (parameters->GetDoPartialVolume()) { parameters->SetXDim(parameters->GetXDim() * RESAMPLING_FACTOR); parameters->SetYDim(parameters->GetYDim() * RESAMPLING_FACTOR); parameters->SetZDim(parameters->GetZDim() * RESAMPLING_FACTOR); parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() * RESAMPLING_FACTOR); parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() / RESAMPLING_FACTOR); } std::mt19937 randomNumberGenerator; std::random_device randomDevice; if (parameters->GetUseRngSeed()) { randomNumberGenerator.seed(parameters->GetRngSeed()); } 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()); } } auto generatedVolume = mitk::pa::InSilicoTissueVolume::New(parameters, &randomNumberGenerator); double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2; if(parameters->GetForceVesselsMoveAlongYDirection()) DIRECTION_VECTOR_INITIAL_VARIANCE = 0; std::uniform_int_distribution randomNumVesselDistribution( parameters->GetMinNumberOfVessels(), parameters->GetMaxNumberOfVessels()); std::uniform_real_distribution randomBendingDistribution( parameters->GetMinVesselBending(), parameters->GetMaxVesselBending()); std::uniform_real_distribution randomAbsorptionDistribution( parameters->GetMinVesselAbsorption(), parameters->GetMaxVesselAbsorption()); std::uniform_real_distribution randomRadiusDistribution( parameters->GetMinVesselRadiusInMillimeters(), parameters->GetMaxVesselRadiusInMillimeters()); std::uniform_real_distribution randomScatteringDistribution( parameters->GetMinVesselScattering(), parameters->GetMaxVesselScattering()); std::uniform_real_distribution randomAnisotropyDistribution( parameters->GetMinVesselAnisotropy(), parameters->GetMaxVesselAnisotropy()); + std::uniform_real_distribution randomOxygenationDistribution(0, 1); std::uniform_int_distribution borderTypeDistribution(0, 3); int numberOfBloodVessels = randomNumVesselDistribution(randomNumberGenerator); generatedVolume->AddIntProperty("numberOfVessels", numberOfBloodVessels); generatedVolume->AddIntProperty("bifurcationFrequency", parameters->GetVesselBifurcationFrequency()); + PropertyCalculator::Pointer m_PropertyCalculator = PropertyCalculator::New(); + MITK_INFO << "Simulating " << numberOfBloodVessels << " vessel structures"; for (int vesselNumber = 0; vesselNumber < numberOfBloodVessels; vesselNumber++) { Vector::Pointer initialPosition = Vector::New(); Vector::Pointer initialDirection = Vector::New(); double initialRadius = randomRadiusDistribution(randomNumberGenerator) / parameters->GetVoxelSpacingInCentimeters() / 10; std::stringstream radiusString; radiusString << "vessel_" << vesselNumber + 1 << "_radius"; generatedVolume->AddDoubleProperty(radiusString.str(), initialRadius); - double absorptionCoefficient = randomAbsorptionDistribution(randomNumberGenerator); - std::stringstream absorptionString; - absorptionString << "vessel_" << vesselNumber + 1 << "_absorption"; - generatedVolume->AddDoubleProperty(absorptionString.str(), absorptionCoefficient); - double bendingFactor = randomBendingDistribution(randomNumberGenerator); std::stringstream bendingString; bendingString << "vessel_" << vesselNumber + 1 << "_bendingFactor"; generatedVolume->AddDoubleProperty(bendingString.str(), bendingFactor); - double vesselScattering = randomScatteringDistribution(randomNumberGenerator); + double absorptionCoefficient = 0; + double vesselScattering = 0; + double vesselAnisotropy = 0; + + if(parameters->GetRandomizeVesselOpticalParams()) + { + double oxygenation = randomOxygenationDistribution(randomNumberGenerator); + int wavelength = parameters->GetTargetWavelength(); + + MITK_INFO << "Randomize vessel parameters. Oxy: " << int(oxygenation*100) << "% @" << wavelength << "nm"; + + auto result = + m_PropertyCalculator->CalculatePropertyForSpecificWavelength(PropertyCalculator::TissueType::BLOOD, wavelength, oxygenation); + + absorptionCoefficient = result.mua; + vesselScattering = result.mus; + vesselAnisotropy = result.g; + + std::stringstream oxygenationString; + oxygenationString << "vessel_" << vesselNumber + 1 << "_oxygenation"; + generatedVolume->AddDoubleProperty(oxygenationString.str(), oxygenation); + } + else + { + absorptionCoefficient = randomAbsorptionDistribution(randomNumberGenerator); + vesselScattering = randomScatteringDistribution(randomNumberGenerator); + vesselAnisotropy = randomAnisotropyDistribution(randomNumberGenerator); + + } + + MITK_INFO << "Setting meta data"; + std::stringstream absorptionString; + absorptionString << "vessel_" << vesselNumber + 1 << "_absorption"; + generatedVolume->AddDoubleProperty(absorptionString.str(), absorptionCoefficient); std::stringstream scatteringString; scatteringString << "vessel_" << vesselNumber + 1 << "_scattering"; generatedVolume->AddDoubleProperty(scatteringString.str(), vesselScattering); - - double vesselAnisotropy = randomAnisotropyDistribution(randomNumberGenerator); std::stringstream anisotropyString; anisotropyString << "vessel_" << vesselNumber + 1 << "_anisotropy"; generatedVolume->AddDoubleProperty(anisotropyString.str(), vesselAnisotropy); + MITK_INFO << "Setting initial direction"; /*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); if(parameters->GetForceVesselsMoveAlongYDirection()) borderType = 2; switch (borderType) { case 0: initialPosition->Randomize(0, 0, initialRadius, parameters->GetYDim() - initialRadius, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &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(parameters->GetXDim() - 1, parameters->GetXDim() - 1, initialRadius, parameters->GetYDim() - initialRadius, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &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, parameters->GetXDim() - initialRadius, 0, 0, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &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, parameters->GetXDim() - initialRadius, parameters->GetYDim() - 1, parameters->GetYDim() - 1, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -2, -1, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; } initialDirection->Normalize(); MITK_INFO << initialPosition->GetElement(0) << " | " << initialPosition->GetElement(1) << " | " << initialPosition->GetElement(2); MITK_INFO << initialDirection->GetElement(0) << " | " << initialDirection->GetElement(1) << " | " << initialDirection->GetElement(2); VesselProperties::Pointer vesselParams = VesselProperties::New(); vesselParams->SetDirectionVector(initialDirection); vesselParams->SetPositionVector(initialPosition); vesselParams->SetRadiusInVoxel(initialRadius); vesselParams->SetAbsorptionCoefficient(absorptionCoefficient); vesselParams->SetScatteringCoefficient(vesselScattering); vesselParams->SetAnisotopyCoefficient(vesselAnisotropy); vesselParams->SetBifurcationFrequency(parameters->GetVesselBifurcationFrequency()); vesselParams->SetDoPartialVolume(parameters->GetDoPartialVolume()); VesselTree::Pointer vesselTree = VesselTree::New(vesselParams); while (!vesselTree->IsFinished()) { vesselTree->Step(generatedVolume, parameters->GetCalculateNewVesselPositionCallback(), bendingFactor, &randomNumberGenerator); } } if (parameters->GetDoPartialVolume()) { VolumeManipulator::RescaleImage(generatedVolume, (1.0 / RESAMPLING_FACTOR)); parameters->SetXDim(parameters->GetXDim() / RESAMPLING_FACTOR); parameters->SetYDim(parameters->GetYDim() / RESAMPLING_FACTOR); parameters->SetZDim(parameters->GetZDim() / RESAMPLING_FACTOR); parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() / RESAMPLING_FACTOR); parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() * RESAMPLING_FACTOR); } generatedVolume->FinalizeVolume(); return generatedVolume; } mitk::pa::InSilicoTissueGenerator::InSilicoTissueGenerator() { } mitk::pa::InSilicoTissueGenerator::~InSilicoTissueGenerator() { } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp index 03922d5835..8ff8e0f1f8 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp @@ -1,155 +1,154 @@ /*=================================================================== 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 "mitkPAPropertyCalculator.h" // us #include #include #include #include #include mitk::pa::PropertyCalculator::Properties mitk::pa::PropertyCalculator:: CalculatePropertyForSpecificWavelength( TissueType tissueType, int wavelength, double bloodOxygenInFraction) { Properties returnValue; if (!m_Valid) { mitkThrow() << "PhotoacousticPropertyGenerator was not loaded properly."; return returnValue; } double bloodOxygenation = bloodOxygenInFraction; double bloodVolumeFraction; double waterVolumeFraction; double fatVolumeFraction; double melanosomesVolumeFraction; double musp500; double fray; double bmie; double g; switch (tissueType) { case AIR: returnValue.mua = 1e-4; returnValue.mus = 1.0; returnValue.g = 1.0; return returnValue; case BLOOD: bloodVolumeFraction = 1.0; waterVolumeFraction = 0.95; fatVolumeFraction = 0.02; melanosomesVolumeFraction = 0; musp500 = 10; fray = 0; bmie = 1; g = 0.9; break; case EPIDERMIS: bloodVolumeFraction = 0; bloodOxygenation = 0; waterVolumeFraction = 0.75; fatVolumeFraction = 0.01; melanosomesVolumeFraction = 0.02; musp500 = 40; fray = 0; bmie = 1; g = 0.9; break; case FAT: bloodVolumeFraction = 0.01; bloodOxygenation = 0.75; waterVolumeFraction = 0.5; fatVolumeFraction = 0.6; melanosomesVolumeFraction = 0; musp500 = 17.47; fray = 0.2; bmie = 0.5; g = 0.9; break; case STANDARD_TISSUE: default: bloodVolumeFraction = 0.02; - bloodOxygenation = 0.75; waterVolumeFraction = 0.8; fatVolumeFraction = 0.05; melanosomesVolumeFraction = 0; musp500 = 25; fray = 0.25; bmie = 1; g = 0.9; break; } // We want the reduced scattering coefficient directly. double musp = musp500 * (fray * pow(wavelength / 500.0, -4.0) + ((1 - fray) * pow(wavelength / 500.0, -bmie))); returnValue.mus = musp; returnValue.mus = 15;//musp; double mua = bloodVolumeFraction*bloodOxygenation*m_SpectralLibMap[MapType::OXYGENATED][wavelength] + bloodVolumeFraction*(1 - bloodOxygenation)*m_SpectralLibMap[MapType::DEOXYGENATED][wavelength] + waterVolumeFraction*m_SpectralLibMap[MapType::WATER][wavelength] + fatVolumeFraction*m_SpectralLibMap[MapType::FATTY][wavelength] + melanosomesVolumeFraction*m_SpectralLibMap[MapType::MELANIN][wavelength]; returnValue.mua = mua; returnValue.g = g; return returnValue; } mitk::pa::PropertyCalculator::PropertyCalculator() { us::ModuleResource resource = us::GetModuleContext()->GetModule()->GetResource("spectralLIB.dat"); if (resource.IsValid()) { us::ModuleResourceStream resourceStream(resource); std::string line; int wavelength = 300; int counter = 0; while (std::getline(resourceStream, line, '\t')) { int currentLineIdx = counter % 6; if (currentLineIdx == 0) wavelength = std::stoi(line); else { std::istringstream lineStream(line); double tempDouble; lineStream >> tempDouble; m_SpectralLibMap[currentLineIdx][wavelength] = tempDouble; } ++counter; } } else { m_Valid = false; } m_Valid = true; } mitk::pa::PropertyCalculator::~PropertyCalculator() { m_SpectralLibMap.clear(); m_Valid = false; } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPATissueGeneratorParameters.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPATissueGeneratorParameters.cpp index bcbeed0f2f..75c3c11e84 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPATissueGeneratorParameters.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPATissueGeneratorParameters.cpp @@ -1,80 +1,83 @@ /*=================================================================== 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 "mitkPATissueGeneratorParameters.h" mitk::pa::TissueGeneratorParameters::TissueGeneratorParameters() { m_XDim = 50; m_YDim = 50; m_ZDim = 50; m_VoxelSpacingInCentimeters = 1; m_DoPartialVolume = false; m_UseRngSeed = false; m_RngSeed = 1337L; m_RandomizePhysicalProperties = false; m_RandomizePhysicalPropertiesPercentage = 0; m_ForceVesselsMoveAlongYDirection = false; m_MinBackgroundAbsorption = 0.1; m_MaxBackgroundAbsorption = 0.1; m_BackgroundScattering = 15; m_BackgroundAnisotropy = 0.9; m_AirAbsorption = 0.0001; m_AirScattering = 1; m_AirAnisotropy = 1; m_AirThicknessInMillimeters = 0; m_SkinAbsorption = 0.1; m_SkinScattering = 15; m_SkinAnisotropy = 0.9; m_SkinThicknessInMillimeters = 0; m_CalculateNewVesselPositionCallback = &VesselMeanderStrategy::CalculateRandomlyDivergingPosition; m_MinNumberOfVessels = 0; m_MaxNumberOfVessels = 0; m_MinVesselBending = 0; m_MaxVesselBending = 0.1; m_MinVesselAbsorption = 1; m_MaxVesselAbsorption = 8; m_MinVesselRadiusInMillimeters = 1; m_MaxVesselRadiusInMillimeters = 3; m_VesselBifurcationFrequency = 25; m_MinVesselScattering = 15; m_MaxVesselScattering = 15; m_MinVesselAnisotropy = 0.9; m_MaxVesselAnisotropy = 0.9; m_MinVesselZOrigin = 10; m_MaxVesselZOrigin = 40; m_MCflag = 1; m_MCLaunchflag = 0; m_MCBoundaryflag = 2; m_MCLaunchPointX = 25; m_MCLaunchPointY = 25; m_MCLaunchPointZ = 2; m_MCFocusPointX = 25; m_MCFocusPointY = 25; m_MCFocusPointZ = 25; m_MCTrajectoryVectorX = 0; m_MCTrajectoryVectorY = 0; m_MCTrajectoryVectorZ = 1; m_MCRadius = 2; m_MCWaist = 4; + + m_RandomizeVesselOpticalParams = false; + m_TargetWavelength = 800; } mitk::pa::TissueGeneratorParameters::~TissueGeneratorParameters() { }