diff --git a/Modules/PhotoacousticsLib/CMakeLists.txt b/Modules/PhotoacousticsLib/CMakeLists.txt index c8a54e6776..0ae41f8f5c 100644 --- a/Modules/PhotoacousticsLib/CMakeLists.txt +++ b/Modules/PhotoacousticsLib/CMakeLists.txt @@ -1,14 +1,15 @@ 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(test) diff --git a/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/CMakeLists.txt b/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/CMakeLists.txt new file mode 100644 index 0000000000..8c8dcb7fba --- /dev/null +++ b/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/CMakeLists.txt @@ -0,0 +1,11 @@ +OPTION(BUILD_PhotoacousticPhantomGenerator "Build MiniApp for generating a PA phantom in silico" ON) + +IF(BUILD_PhotoacousticPhantomGenerator) + PROJECT( MitkPAPhantomGenerator ) + mitk_create_executable(PAPhantomGenerator + DEPENDS MitkCommandLine MitkCore MitkPhotoacousticsLib + PACKAGE_DEPENDS + CPP_FILES PAPhantomGenerator.cpp) + + install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) + ENDIF() diff --git a/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/PAPhantomGenerator.cpp b/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/PAPhantomGenerator.cpp new file mode 100644 index 0000000000..decd6880be --- /dev/null +++ b/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/PAPhantomGenerator.cpp @@ -0,0 +1,216 @@ +/*=================================================================== + +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 + +using namespace mitk::pa; + +TissueGeneratorParameters::Pointer CreatePhantom_08_03_18_Parameters() +{ + auto returnParameters = TissueGeneratorParameters::New(); + returnParameters->SetAirThicknessInMillimeters(12); + returnParameters->SetMinBackgroundAbsorption(0.1); + returnParameters->SetMaxBackgroundAbsorption(0.1); + returnParameters->SetBackgroundAnisotropy(0.9); + returnParameters->SetBackgroundScattering(15); + returnParameters->SetCalculateNewVesselPositionCallback(&VesselMeanderStrategy::CalculateNewPositionInStraightLine); + returnParameters->SetDoPartialVolume(false); + returnParameters->SetMinNumberOfVessels(5); + returnParameters->SetMaxNumberOfVessels(5); + returnParameters->SetMinVesselAbsorption(10); + returnParameters->SetMaxVesselAbsorption(10); + returnParameters->SetMinVesselAnisotropy(0.9); + returnParameters->SetMaxVesselAnisotropy(0.9); + returnParameters->SetMinVesselBending(0.1); + returnParameters->SetMaxVesselBending(0.3); + returnParameters->SetMinVesselRadiusInMillimeters(0.5); + returnParameters->SetMaxVesselRadiusInMillimeters(4); + returnParameters->SetMinVesselScattering(15); + returnParameters->SetMaxVesselScattering(15); + returnParameters->SetMinVesselZOrigin(2.2); + returnParameters->SetMaxVesselZOrigin(4); + returnParameters->SetVesselBifurcationFrequency(5000); + returnParameters->SetRandomizePhysicalProperties(false); + returnParameters->SetSkinThicknessInMillimeters(0); + returnParameters->SetUseRngSeed(false); + returnParameters->SetVoxelSpacingInCentimeters(0.0075); + returnParameters->SetXDim(560); + returnParameters->SetYDim(800); + returnParameters->SetZDim(720); + 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 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; +} + +int main(int argc, char * argv[]) +{ + auto input = parseInput(argc, argv); + double absorptions[10] = {1, 3, 5, 7, 9, 11, 13, 15, 17, 19}; + unsigned int iterationNumber = 0; + + for(;iterationNumber<10; iterationNumber++) + { + + auto parameters = CreatePhantom_08_03_18_Parameters(); + parameters->SetMinVesselAbsorption(absorptions[iterationNumber]); + parameters->SetMaxVesselAbsorption(absorptions[iterationNumber]); + MITK_INFO(input.verbose) << "Generating tissue.."; + auto resultTissue = PhantomTissueGenerator::GeneratePhantomData(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/Phantom_" + input.identifyer + + "_" + std::to_string(iterationNumber) + ".nrrd"; + mitk::IOUtil::Save(resultTissue->ConvertToMitkImage(), savePath); + std::string outputPath = input.saveFolderPath + "output/Phantom_" + input.identifyer + + "_" + std::to_string(iterationNumber) + "/"; + + if (!itksys::SystemTools::FileIsDirectory(outputPath)) + { + itksys::SystemTools::MakeDirectory(outputPath); + } + + outputPath = outputPath + "Fluence_Phantom_" + input.identifyer + "_" + std::to_string(iterationNumber); + + MITK_INFO(input.verbose) << "Simulating fluence.."; + + int result = -4; + if(!input.probePath.empty()) + result = std::system(std::string(input.exePath + " -i " + savePath + " -o " + + (outputPath + ".nrrd") + + " -yo " + "0" + " -p " + input.probePath + + " -n 50000").c_str()); + else + result = std::system(std::string(input.exePath + " -i " + savePath + " -o " + + (outputPath + ".nrrd") + + " -yo " + "0" + " -n 50000").c_str()); + + MITK_INFO << result; + MITK_INFO(input.verbose) << "Simulating fluence..[Done]"; + } + +} diff --git a/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/files.cmake b/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/files.cmake new file mode 100644 index 0000000000..07f708562d --- /dev/null +++ b/Modules/PhotoacousticsLib/MitkPAPhantomGenerator/files.cmake @@ -0,0 +1,3 @@ +set(CPP_FILES + PAPhantomGenerator.cpp +) diff --git a/Modules/PhotoacousticsLib/files.cmake b/Modules/PhotoacousticsLib/files.cmake index 28e6b805ad..3f62323168 100644 --- a/Modules/PhotoacousticsLib/files.cmake +++ b/Modules/PhotoacousticsLib/files.cmake @@ -1,54 +1,56 @@ SET(H_FILES include/mitkPAPropertyCalculator.h include/mitkPAVector.h include/mitkPATissueGeneratorParameters.h include/mitkPAInSilicoTissueVolume.h + include/mitkPAPhantomTissueGenerator.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/mitkPAPhantomTissueGenerator.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/mitkPAPhantomTissueGenerator.h b/Modules/PhotoacousticsLib/include/mitkPAPhantomTissueGenerator.h new file mode 100644 index 0000000000..f60c62de7c --- /dev/null +++ b/Modules/PhotoacousticsLib/include/mitkPAPhantomTissueGenerator.h @@ -0,0 +1,52 @@ +/*=================================================================== + +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 mitkPhotoacousticPhantomTissueGenerator_h +#define mitkPhotoacousticPhantomTissueGenerator_h + +#include +#include +#include +#include + +#include +#include + +#include "mitkPAVesselTree.h" +#include "mitkPAInSilicoTissueVolume.h" + +#include "mitkCommon.h" + +namespace mitk { + namespace pa { + class MITKPHOTOACOUSTICSLIB_EXPORT PhantomTissueGenerator final + { + public: + + /** + * @brief GenerateInSilicoData This method will return a InSilicoTissueVolume created in terms of the given parameters. + * @param parameters + * @return + */ + static InSilicoTissueVolume::Pointer GeneratePhantomData(TissueGeneratorParameters::Pointer parameters); + + private: + PhantomTissueGenerator(); + virtual ~PhantomTissueGenerator(); + }; + } +} +#endif diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPAPhantomTissueGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPAPhantomTissueGenerator.cpp new file mode 100644 index 0000000000..5b264c08ae --- /dev/null +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPAPhantomTissueGenerator.cpp @@ -0,0 +1,158 @@ +/*=================================================================== + +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 "mitkPAPhantomTissueGenerator.h" +#include "mitkPAVector.h" +#include "mitkPAVolumeManipulator.h" + +mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::PhantomTissueGenerator::GeneratePhantomData( + 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); + } + + parameters->SetVesselBifurcationFrequency(10000); + + 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); + + const unsigned int NUMER_OF_VESSELS = 5; + const double START_DEPTH_IN_CM = 2.10; + const double START_X_IN_CM = 0.76; + const double RADIUS_IN_MM = 0.6125; + const double INCREMENT_XZ_IN_CM = 0.50; + double ABSORPTION_PER_CM = parameters->GetMinVesselAbsorption(); + + generatedVolume->AddIntProperty("numberOfVessels", NUMER_OF_VESSELS); + generatedVolume->AddIntProperty("bifurcationFrequency", parameters->GetVesselBifurcationFrequency()); + + for (unsigned int vesselNumber = 0; vesselNumber < NUMER_OF_VESSELS; vesselNumber++) + { + Vector::Pointer initialPosition = Vector::New(); + Vector::Pointer initialDirection = Vector::New(); + + double initialRadius = RADIUS_IN_MM / parameters->GetVoxelSpacingInCentimeters() / 10; + std::stringstream radiusString; + radiusString << "vessel_" << vesselNumber + 1 << "_radius"; + generatedVolume->AddDoubleProperty(radiusString.str(), initialRadius); + + double absorptionCoefficient = ABSORPTION_PER_CM; + std::stringstream absorptionString; + absorptionString << "vessel_" << vesselNumber + 1 << "_absorption"; + generatedVolume->AddDoubleProperty(absorptionString.str(), absorptionCoefficient); + + double bendingFactor = 0; + std::stringstream bendingString; + bendingString << "vessel_" << vesselNumber + 1 << "_bendingFactor"; + generatedVolume->AddDoubleProperty(bendingString.str(), bendingFactor); + + double vesselScattering = 15; + std::stringstream scatteringString; + scatteringString << "vessel_" << vesselNumber + 1 << "_scattering"; + generatedVolume->AddDoubleProperty(scatteringString.str(), vesselScattering); + + double vesselAnisotropy = 0.9; + std::stringstream anisotropyString; + anisotropyString << "vessel_" << vesselNumber + 1 << "_anisotropy"; + generatedVolume->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) + */ + + double zposition = (START_DEPTH_IN_CM + (vesselNumber*INCREMENT_XZ_IN_CM)) / parameters->GetVoxelSpacingInCentimeters(); + + double xposition = (START_X_IN_CM + (vesselNumber*INCREMENT_XZ_IN_CM)) / parameters->GetVoxelSpacingInCentimeters(); + + + initialPosition->Randomize(xposition, xposition, 0, 0, zposition, zposition, &randomNumberGenerator); + initialDirection->Randomize(0, 0, 1, 1, 0, 0, &randomNumberGenerator); + + 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::PhantomTissueGenerator::PhantomTissueGenerator() +{ +} + +mitk::pa::PhantomTissueGenerator::~PhantomTissueGenerator() +{ +}