diff --git a/Modules/PhotoacousticsLib/MitkTissueBatchGenerator/TissueBatchGenerator.cpp b/Modules/PhotoacousticsLib/MitkTissueBatchGenerator/TissueBatchGenerator.cpp index 462e872e2b..61d6b94ab9 100644 --- a/Modules/PhotoacousticsLib/MitkTissueBatchGenerator/TissueBatchGenerator.cpp +++ b/Modules/PhotoacousticsLib/MitkTissueBatchGenerator/TissueBatchGenerator.cpp @@ -1,224 +1,394 @@ /*=================================================================== 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 CreateMultiHB_13_02_18_Parameters() +{ + auto returnParameters = TissueGeneratorParameters::New(); + returnParameters->SetAirThicknessInMillimeters(1.8); + returnParameters->SetMinBackgroundAbsorption(0.001); + returnParameters->SetMaxBackgroundAbsorption(0.2); + returnParameters->SetBackgroundAnisotropy(0.9); + returnParameters->SetBackgroundScattering(15); + returnParameters->SetCalculateNewVesselPositionCallback(&VesselMeanderStrategy::CalculateRandomlyDivergingPosition); + returnParameters->SetDoPartialVolume(true); + returnParameters->SetMinNumberOfVessels(1); + returnParameters->SetMaxNumberOfVessels(7); + returnParameters->SetMinVesselAbsorption(1); + returnParameters->SetMaxVesselAbsorption(12); + returnParameters->SetMinVesselAnisotropy(0.9); + returnParameters->SetMaxVesselAnisotropy(0.9); + returnParameters->SetMinVesselBending(0); + returnParameters->SetMaxVesselBending(0.2); + returnParameters->SetMinVesselRadiusInMillimeters(0.5); + returnParameters->SetMaxVesselRadiusInMillimeters(6); + returnParameters->SetMinVesselScattering(15); + returnParameters->SetMaxVesselScattering(15); + returnParameters->SetMinVesselZOrigin(1); + returnParameters->SetMaxVesselZOrigin(3); + returnParameters->SetVesselBifurcationFrequency(5000); + returnParameters->SetRandomizePhysicalProperties(false); + returnParameters->SetSkinThicknessInMillimeters(0); + returnParameters->SetUseRngSeed(false); + returnParameters->SetVoxelSpacingInCentimeters(0.06); + returnParameters->SetXDim(70); + returnParameters->SetYDim(100); + returnParameters->SetZDim(100); + returnParameters->SetMCflag(4); + return returnParameters; +} + +TissueGeneratorParameters::Pointer CreateBaselineHB_13_02_18_Parameters() +{ + auto returnParameters = TissueGeneratorParameters::New(); + returnParameters->SetAirThicknessInMillimeters(1.8); + returnParameters->SetMinBackgroundAbsorption(0.001); + returnParameters->SetMaxBackgroundAbsorption(0.2); + returnParameters->SetBackgroundAnisotropy(0.9); + returnParameters->SetBackgroundScattering(15); + returnParameters->SetCalculateNewVesselPositionCallback(&VesselMeanderStrategy::CalculateRandomlyDivergingPosition); + returnParameters->SetDoPartialVolume(true); + returnParameters->SetMinNumberOfVessels(1); + returnParameters->SetMaxNumberOfVessels(1); + returnParameters->SetMinVesselAbsorption(4.73); + returnParameters->SetMaxVesselAbsorption(4.73); + returnParameters->SetMinVesselAnisotropy(0.9); + returnParameters->SetMaxVesselAnisotropy(0.9); + returnParameters->SetMinVesselBending(0); + returnParameters->SetMaxVesselBending(0.2); + returnParameters->SetMinVesselRadiusInMillimeters(3); + returnParameters->SetMaxVesselRadiusInMillimeters(3); + returnParameters->SetMinVesselScattering(15); + returnParameters->SetMaxVesselScattering(15); + returnParameters->SetMinVesselZOrigin(1); + returnParameters->SetMaxVesselZOrigin(3); + returnParameters->SetVesselBifurcationFrequency(5000); + returnParameters->SetRandomizePhysicalProperties(false); + returnParameters->SetSkinThicknessInMillimeters(0); + returnParameters->SetUseRngSeed(false); + returnParameters->SetVoxelSpacingInCentimeters(0.06); + returnParameters->SetXDim(70); + returnParameters->SetYDim(100); + returnParameters->SetZDim(100); + returnParameters->SetMCflag(4); + return returnParameters; +} + +TissueGeneratorParameters::Pointer CreateSingleVesselHeterogeneousBackground_08_02_18_Parameters() +{ + auto returnParameters = TissueGeneratorParameters::New(); + returnParameters->SetAirThicknessInMillimeters(1.8); + returnParameters->SetMinBackgroundAbsorption(0.001); + returnParameters->SetMaxBackgroundAbsorption(0.2); + returnParameters->SetBackgroundAnisotropy(0.9); + returnParameters->SetBackgroundScattering(15); + returnParameters->SetCalculateNewVesselPositionCallback(&VesselMeanderStrategy::CalculateRandomlyDivergingPosition); + returnParameters->SetDoPartialVolume(true); + returnParameters->SetMinNumberOfVessels(1); + returnParameters->SetMaxNumberOfVessels(1); + returnParameters->SetMinVesselAbsorption(1); + returnParameters->SetMaxVesselAbsorption(12); + returnParameters->SetMinVesselAnisotropy(0.9); + returnParameters->SetMaxVesselAnisotropy(0.9); + returnParameters->SetMinVesselBending(0); + returnParameters->SetMaxVesselBending(0.2); + returnParameters->SetMinVesselRadiusInMillimeters(0.5); + returnParameters->SetMaxVesselRadiusInMillimeters(6); + returnParameters->SetMinVesselScattering(15); + returnParameters->SetMaxVesselScattering(15); + returnParameters->SetMinVesselZOrigin(1); + returnParameters->SetMaxVesselZOrigin(3); + returnParameters->SetVesselBifurcationFrequency(5000); + returnParameters->SetRandomizePhysicalProperties(false); + returnParameters->SetSkinThicknessInMillimeters(0); + returnParameters->SetUseRngSeed(false); + returnParameters->SetVoxelSpacingInCentimeters(0.06); + returnParameters->SetXDim(70); + returnParameters->SetYDim(100); + returnParameters->SetZDim(100); + returnParameters->SetMCflag(4); + return returnParameters; +} + +TissueGeneratorParameters::Pointer CreateMultivessel_19_12_17_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::CalculateRandomlyDivergingPosition); + returnParameters->SetDoPartialVolume(true); + returnParameters->SetMinNumberOfVessels(1); + returnParameters->SetMaxNumberOfVessels(7); + returnParameters->SetMinVesselAbsorption(2); + returnParameters->SetMaxVesselAbsorption(8); + 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.06); + returnParameters->SetXDim(70); + returnParameters->SetYDim(100); + returnParameters->SetZDim(100); + return returnParameters; +} + TissueGeneratorParameters::Pointer CreateMultivessel_19_10_17_Parameters() { auto returnParameters = TissueGeneratorParameters::New(); returnParameters->SetAirThicknessInMillimeters(12); - returnParameters->SetBackgroundAbsorption(0.1); + returnParameters->SetMinBackgroundAbsorption(0.1); + returnParameters->SetMaxBackgroundAbsorption(0.1); returnParameters->SetBackgroundAnisotropy(0.9); returnParameters->SetBackgroundScattering(15); returnParameters->SetCalculateNewVesselPositionCallback(&VesselMeanderStrategy::CalculateRandomlyDivergingPosition); returnParameters->SetDoPartialVolume(true); returnParameters->SetMinNumberOfVessels(1); returnParameters->SetMaxNumberOfVessels(7); returnParameters->SetMinVesselAbsorption(2); returnParameters->SetMaxVesselAbsorption(8); 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.03); returnParameters->SetXDim(140); returnParameters->SetYDim(200); returnParameters->SetZDim(180); return returnParameters; } TissueGeneratorParameters::Pointer CreateSinglevessel_19_10_17_Parameters() { auto returnParameters = TissueGeneratorParameters::New(); returnParameters->SetAirThicknessInMillimeters(12); - returnParameters->SetBackgroundAbsorption(0.1); + returnParameters->SetMinBackgroundAbsorption(0.1); + returnParameters->SetMaxBackgroundAbsorption(0.1); returnParameters->SetBackgroundAnisotropy(0.9); returnParameters->SetBackgroundScattering(15); returnParameters->SetCalculateNewVesselPositionCallback(&VesselMeanderStrategy::CalculateRandomlyDivergingPosition); returnParameters->SetDoPartialVolume(true); returnParameters->SetMinNumberOfVessels(1); returnParameters->SetMaxNumberOfVessels(1); returnParameters->SetMinVesselAbsorption(2); returnParameters->SetMaxVesselAbsorption(8); 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.03); returnParameters->SetXDim(140); returnParameters->SetYDim(200); returnParameters->SetZDim(180); 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(), false); - parser.endGroup(); - parser.beginGroup("Optional parameters"); + 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); unsigned int iterationNumber = 0; while (true) { - auto parameters = CreateSinglevessel_19_10_17_Parameters(); + auto parameters = CreateBaselineHB_13_02_18_Parameters(); 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/MultiVessel_" + input.identifyer + "_" + std::to_string(iterationNumber) + ".nrrd"; + std::string savePath = input.saveFolderPath + "input/BaselineHB_" + input.identifyer + + "_" + std::to_string(iterationNumber) + ".nrrd"; mitk::IOUtil::Save(resultTissue->ConvertToMitkImage(), savePath); - std::string outputPath = input.saveFolderPath + "output/Fluence_MultiVessel_" + input.identifyer + "_" + std::to_string(iterationNumber) + ".nrrd"; + std::string outputPath = input.saveFolderPath + "output/BaselineHB_" + input.identifyer + + "_" + std::to_string(iterationNumber) + "/"; + + if (!itksys::SystemTools::FileIsDirectory(outputPath)) + { + itksys::SystemTools::MakeDirectory(outputPath); + } + + outputPath = outputPath + "Fluence_BaselineHB_" + input.identifyer + "_" + std::to_string(iterationNumber); MITK_INFO(input.verbose) << "Simulating fluence.."; - int result = std::system(std::string(input.exePath + " -i " + savePath + " -o " + outputPath + " -p " + input.probePath + " -n 100000000").c_str()); - MITK_INFO(input.verbose) << "Simulating fluence"<< result <<"..[Done]"; + for(double yo = -1.8; yo <= 1.81; yo=yo+0.12) + { + 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 + "_yo" + yo_string + ".nrrd") + + " -yo " + yo_string + " -p " + input.probePath + + " -n 100000000").c_str()); + else + result = std::system(std::string(input.exePath + " -i " + savePath + " -o " + + (outputPath + "_yo" + yo_string + ".nrrd") + + " -yo " + yo_string + " -n 100000000").c_str()); + MITK_INFO << "yo: " << yo_string << ": " << result; + } + + MITK_INFO(input.verbose) << "Simulating fluence..[Done]"; + iterationNumber++; } } diff --git a/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h b/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h index 2b3575fbb1..de81cb5520 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h +++ b/Modules/PhotoacousticsLib/include/mitkPAInSilicoTissueVolume.h @@ -1,160 +1,164 @@ /*=================================================================== 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 MITKPHOTOACOUSTICVOLUME_H #define MITKPHOTOACOUSTICVOLUME_H #include #include #include #include #include //Includes for smart pointer usage #include "mitkCommon.h" #include "itkLightObject.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT InSilicoTissueVolume : public itk::LightObject { public: mitkClassMacroItkParent(InSilicoTissueVolume, itk::LightObject) - mitkNewMacro1Param(Self, TissueGeneratorParameters::Pointer) + mitkNewMacro2Param(Self, TissueGeneratorParameters::Pointer, std::mt19937*) - enum SegmentationType + enum SegmentationType { AIR = -1, BACKGROUND = 0, VESSEL = 1, FAT = 2, SKIN = 3 }; /** * @brief ConvertToMitkImage * @return a pointer to an mitk image containing this volume. */ mitk::Image::Pointer ConvertToMitkImage(); /** * @brief SetVolumeValues sets the values for aborption, scattering and anisotropy at the specified voxel location. * * @param x * @param y * @param z * @param absorption * @param scattering * @param anisotropy * @param segmentType */ void SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy, SegmentationType segmentType); /** * @brief SetVolumeValues sets the values for aborption, scattering and anisotropy at the specified voxel location. * * @param x * @param y * @param z * @param absorption * @param scattering * @param anisotropy */ void SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy); /** * @brief IsInsideVolume * * @param x * @param y * @param z * @return true if the voxel location is inside the volume */ bool IsInsideVolume(int x, int y, int z); /** * @brief AddDoubleProperty adds a persistent property to the volume, which will be exported to the mitk image. * * @param label * @param value */ void AddDoubleProperty(std::string label, double value); /** * @brief AddIntProperty adds a persistent property to the volume, which will be exported to the mitk image. * * @param label * @param value */ void AddIntProperty(std::string label, int value); Volume::Pointer GetAbsorptionVolume(); Volume::Pointer GetScatteringVolume(); Volume::Pointer GetAnisotropyVolume(); Volume::Pointer GetSegmentationVolume(); void SetAbsorptionVolume(Volume::Pointer volume); void SetScatteringVolume(Volume::Pointer volume); void SetAnisotropyVolume(Volume::Pointer volume); void SetSegmentationVolume(Volume::Pointer volume); double GetSpacing(); void SetSpacing(double spacing); void FinalizeVolume(); itkGetMacro(TissueParameters, TissueGeneratorParameters::Pointer); itkGetMacro(TDim, unsigned int); static InSilicoTissueVolume::Pointer New(mitk::pa::Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList); protected: - InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters); + InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters, std::mt19937* rng); InSilicoTissueVolume(Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList); virtual ~InSilicoTissueVolume(); mitk::pa::Volume::Pointer m_AbsorptionVolume; mitk::pa::Volume::Pointer m_ScatteringVolume; mitk::pa::Volume::Pointer m_AnisotropyVolume; mitk::pa::Volume::Pointer m_SegmentationVolume; TissueGeneratorParameters::Pointer m_TissueParameters; unsigned int m_TDim; + double m_InitialBackgroundAbsorption; + + std::mt19937* m_Rng; + void RandomizeTissueCoefficients(long rngSeed, bool useRngSeed, double percentage); mitk::PropertyList::Pointer m_PropertyList; private: void FillZLayer(int x, int y, double startIdx, double endIdx, double absorption, double scattering, double anisotropy, SegmentationType segmentationType); void AddSkinAndAirLayers(); void UpdatePropertyList(); }; } } #endif // MITKPHOTOACOUSTICVOLUME_H diff --git a/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h b/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h index 6be09cc42f..946c912eb0 100644 --- a/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h +++ b/Modules/PhotoacousticsLib/include/mitkPATissueGeneratorParameters.h @@ -1,211 +1,214 @@ /*=================================================================== 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(BackgroundAbsorption, double) + 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) 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(BackgroundAbsorption, double) + 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) protected: TissueGeneratorParameters(); virtual ~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; - double m_BackgroundAbsorption; + 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; 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/Domain/Volume/mitkPAInSilicoTissueVolume.cpp b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp index eb8e036a2b..301a5e91f1 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp @@ -1,360 +1,388 @@ /*=================================================================== 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 -mitk::pa::InSilicoTissueVolume::InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters) +mitk::pa::InSilicoTissueVolume::InSilicoTissueVolume(TissueGeneratorParameters::Pointer parameters, std::mt19937* rng) { { unsigned int xDim = parameters->GetXDim(); unsigned int yDim = parameters->GetYDim(); unsigned int zDim = parameters->GetZDim(); m_TDim = 4; unsigned int size = xDim * yDim * zDim; double* absorptionArray = new double[size]; double* scatteringArray = new double[size]; double* anisotropyArray = new double[size]; double* segmentationArray = new double[size]; + m_InitialBackgroundAbsorption = (parameters->GetMinBackgroundAbsorption() + parameters->GetMaxBackgroundAbsorption()) / 2; + m_Rng = rng; + for (unsigned int index = 0; index < size; index++) { - absorptionArray[index] = parameters->GetBackgroundAbsorption(); + absorptionArray[index] = m_InitialBackgroundAbsorption; scatteringArray[index] = parameters->GetBackgroundScattering(); anisotropyArray[index] = parameters->GetBackgroundAnisotropy(); segmentationArray[index] = SegmentationType::BACKGROUND; } m_AbsorptionVolume = Volume::New(absorptionArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); m_ScatteringVolume = Volume::New(scatteringArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); m_AnisotropyVolume = Volume::New(anisotropyArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); m_SegmentationVolume = Volume::New(segmentationArray, xDim, yDim, zDim, parameters->GetVoxelSpacingInCentimeters()); } m_TissueParameters = parameters; m_PropertyList = mitk::PropertyList::New(); UpdatePropertyList(); } void mitk::pa::InSilicoTissueVolume::UpdatePropertyList() { //Set properties AddIntProperty("mcflag", m_TissueParameters->GetMCflag()); AddIntProperty("launchflag", m_TissueParameters->GetMCLaunchflag()); AddIntProperty("boundaryflag", m_TissueParameters->GetMCBoundaryflag()); AddDoubleProperty("launchPointX", m_TissueParameters->GetMCLaunchPointX()); AddDoubleProperty("launchPointY", m_TissueParameters->GetMCLaunchPointY()); AddDoubleProperty("launchPointZ", m_TissueParameters->GetMCLaunchPointZ()); AddDoubleProperty("focusPointX", m_TissueParameters->GetMCFocusPointX()); AddDoubleProperty("focusPointY", m_TissueParameters->GetMCFocusPointY()); AddDoubleProperty("focusPointZ", m_TissueParameters->GetMCFocusPointZ()); AddDoubleProperty("trajectoryVectorX", m_TissueParameters->GetMCTrajectoryVectorX()); AddDoubleProperty("trajectoryVectorY", m_TissueParameters->GetMCTrajectoryVectorY()); AddDoubleProperty("trajectoryVectorZ", m_TissueParameters->GetMCTrajectoryVectorZ()); AddDoubleProperty("radius", m_TissueParameters->GetMCRadius()); AddDoubleProperty("waist", m_TissueParameters->GetMCWaist()); AddDoubleProperty("partialVolume", m_TissueParameters->GetDoPartialVolume()); - AddDoubleProperty("standardTissueAbsorption", m_TissueParameters->GetBackgroundAbsorption()); + AddDoubleProperty("standardTissueAbsorptionMin", m_TissueParameters->GetMinBackgroundAbsorption()); + AddDoubleProperty("standardTissueAbsorptionMax", m_TissueParameters->GetMaxBackgroundAbsorption()); AddDoubleProperty("standardTissueScattering", m_TissueParameters->GetBackgroundScattering()); AddDoubleProperty("standardTissueAnisotropy", m_TissueParameters->GetBackgroundAnisotropy()); AddDoubleProperty("airThickness", m_TissueParameters->GetAirThicknessInMillimeters()); AddDoubleProperty("skinThickness", m_TissueParameters->GetSkinThicknessInMillimeters()); } mitk::pa::InSilicoTissueVolume::InSilicoTissueVolume( Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList) { m_AbsorptionVolume = absorptionVolume; m_ScatteringVolume = scatteringVolume; m_AnisotropyVolume = anisotropyVolume; m_SegmentationVolume = segmentationVolume; m_TissueParameters = tissueParameters; m_PropertyList = propertyList; if (m_SegmentationVolume.IsNotNull()) m_TDim = 4; else m_TDim = 3; } double mitk::pa::InSilicoTissueVolume::GetSpacing() { return m_AbsorptionVolume->GetSpacing(); } void mitk::pa::InSilicoTissueVolume::SetSpacing(double spacing) { m_AbsorptionVolume->SetSpacing(spacing); m_ScatteringVolume->SetSpacing(spacing); m_AnisotropyVolume->SetSpacing(spacing); m_SegmentationVolume->SetSpacing(spacing); } void mitk::pa::InSilicoTissueVolume::AddDoubleProperty(std::string label, double value) { m_PropertyList->SetDoubleProperty(label.c_str(), value); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New(label)); } void mitk::pa::InSilicoTissueVolume::AddIntProperty(std::string label, int value) { m_PropertyList->SetIntProperty(label.c_str(), value); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New(label)); } mitk::Image::Pointer mitk::pa::InSilicoTissueVolume::ConvertToMitkImage() { mitk::Image::Pointer resultImage = mitk::Image::New(); mitk::PixelType TPixel = mitk::MakeScalarPixelType(); unsigned int* dimensionsOfImage = new unsigned int[4]; // Copy dimensions dimensionsOfImage[0] = m_TissueParameters->GetYDim(); dimensionsOfImage[1] = m_TissueParameters->GetXDim(); dimensionsOfImage[2] = m_TissueParameters->GetZDim(); dimensionsOfImage[3] = m_TDim; resultImage->Initialize(TPixel, 4, dimensionsOfImage, 1); mitk::Vector3D spacing; spacing.Fill(m_TissueParameters->GetVoxelSpacingInCentimeters()); resultImage->SetSpacing(spacing); resultImage->SetImportVolume(m_AbsorptionVolume->GetData(), 0, 0, mitk::Image::CopyMemory); resultImage->SetImportVolume(m_ScatteringVolume->GetData(), 1, 0, mitk::Image::CopyMemory); resultImage->SetImportVolume(m_AnisotropyVolume->GetData(), 2, 0, mitk::Image::CopyMemory); resultImage->SetImportVolume(m_SegmentationVolume->GetData(), 3, 0, mitk::Image::CopyMemory); resultImage->SetPropertyList(m_PropertyList); return resultImage; } mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::InSilicoTissueVolume::New( Volume::Pointer absorptionVolume, Volume::Pointer scatteringVolume, Volume::Pointer anisotropyVolume, Volume::Pointer segmentationVolume, TissueGeneratorParameters::Pointer tissueParameters, mitk::PropertyList::Pointer propertyList) { InSilicoTissueVolume::Pointer smartPtr = new InSilicoTissueVolume( absorptionVolume, scatteringVolume, anisotropyVolume, segmentationVolume, tissueParameters, propertyList); smartPtr->UnRegister(); return smartPtr; } mitk::pa::InSilicoTissueVolume::~InSilicoTissueVolume() { m_AbsorptionVolume = nullptr; m_ScatteringVolume = nullptr; m_AnisotropyVolume = nullptr; m_SegmentationVolume = nullptr; m_TissueParameters = nullptr; m_PropertyList = nullptr; } void mitk::pa::InSilicoTissueVolume::SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy) { if (IsInsideVolume(x, y, z)) { m_AbsorptionVolume->SetData(absorption, x, y, z); m_ScatteringVolume->SetData(scattering, x, y, z); m_AnisotropyVolume->SetData(anisotropy, x, y, z); } } void mitk::pa::InSilicoTissueVolume::SetVolumeValues(int x, int y, int z, double absorption, double scattering, double anisotropy, SegmentationType segmentType) { if (IsInsideVolume(x, y, z)) { m_AbsorptionVolume->SetData(absorption, x, y, z); m_ScatteringVolume->SetData(scattering, x, y, z); m_AnisotropyVolume->SetData(anisotropy, x, y, z); m_SegmentationVolume->SetData(segmentType, x, y, z); } } bool mitk::pa::InSilicoTissueVolume::IsInsideVolume(int x, int y, int z) { return x >= 0 && x < m_TissueParameters->GetXDim() && y >= 0 && y < m_TissueParameters->GetYDim() && z >= 0 && z < m_TissueParameters->GetZDim(); } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetAbsorptionVolume() { return m_AbsorptionVolume; } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetSegmentationVolume() { return m_SegmentationVolume; } void mitk::pa::InSilicoTissueVolume::FinalizeVolume() { AddSkinAndAirLayers(); // If specified, randomize all tissue parameters if (m_TissueParameters->GetRandomizePhysicalProperties()) + { RandomizeTissueCoefficients(m_TissueParameters->GetUseRngSeed(), m_TissueParameters->GetRngSeed(), m_TissueParameters->GetRandomizePhysicalPropertiesPercentage()); + } + + unsigned int xDim = m_TissueParameters->GetXDim(); + unsigned int yDim = m_TissueParameters->GetYDim(); + unsigned int zDim = m_TissueParameters->GetZDim(); + + std::uniform_real_distribution randomBackgroundAbsorptionDistribution( + m_TissueParameters->GetMinBackgroundAbsorption(), m_TissueParameters->GetMaxBackgroundAbsorption()); + + + for (unsigned int z = 0; z < zDim; z++) + { + for (unsigned int y = 0; y < yDim; y++) + { + for (unsigned int x = 0; x < xDim; x++) + { + if(fabs(m_AbsorptionVolume->GetData(x, y, z)-m_InitialBackgroundAbsorption) < mitk::eps) + { + m_AbsorptionVolume->SetData(randomBackgroundAbsorptionDistribution(*m_Rng), x, y, z); + } + } + } + } } void mitk::pa::InSilicoTissueVolume::AddSkinAndAirLayers() { //Calculate the index location according to thickness in cm double airvoxel = (m_TissueParameters->GetAirThicknessInMillimeters() / m_TissueParameters->GetVoxelSpacingInCentimeters()) / 10; double skinvoxel = airvoxel + (m_TissueParameters->GetSkinThicknessInMillimeters() / m_TissueParameters->GetVoxelSpacingInCentimeters()) / 10; for (int y = 0; y < m_TissueParameters->GetYDim(); y++) { for (int x = 0; x < m_TissueParameters->GetXDim(); x++) { // Add air from index 0 to airvoxel if (m_TissueParameters->GetAirThicknessInMillimeters() > mitk::eps) { FillZLayer(x, y, 0, airvoxel, m_TissueParameters->GetAirAbsorption(), m_TissueParameters->GetAirScattering(), m_TissueParameters->GetAirAnisotropy(), SegmentationType::AIR); } //Add skin from index airvoxel to skinvoxel if (m_TissueParameters->GetSkinThicknessInMillimeters() > mitk::eps) { FillZLayer(x, y, airvoxel, skinvoxel, m_TissueParameters->GetSkinAbsorption(), m_TissueParameters->GetSkinScattering(), m_TissueParameters->GetSkinAnisotropy(), SegmentationType::SKIN); } } } } void mitk::pa::InSilicoTissueVolume::FillZLayer(int x, int y, double startIdx, double endIdx, double absorption, double scattering, double anisotropy, SegmentationType segmentationType) { for (int z = startIdx; z < endIdx; z++) { if (IsInsideVolume(x, y, z)) { if (endIdx - z < 1) { //Simulate partial volume effects m_AbsorptionVolume->SetData((1 - (endIdx - z)) * m_AbsorptionVolume->GetData(x, y, z) + (endIdx - z) * absorption, x, y, z); m_ScatteringVolume->SetData((1 - (endIdx - z)) * m_ScatteringVolume->GetData(x, y, z) + (endIdx - z) * scattering, x, y, z); m_AnisotropyVolume->SetData((1 - (endIdx - z)) * m_AnisotropyVolume->GetData(x, y, z) + (endIdx - z) * anisotropy, x, y, z); if (endIdx - z > 0.5) { //Only put the segmentation label if more than half of the partial volume is the wanted tissue type m_SegmentationVolume->SetData(segmentationType, x, y, z); } } else { m_AbsorptionVolume->SetData(absorption, x, y, z); m_ScatteringVolume->SetData(scattering, x, y, z); m_AnisotropyVolume->SetData(anisotropy, x, y, z); m_SegmentationVolume->SetData(segmentationType, x, y, z); } } } } void mitk::pa::InSilicoTissueVolume::RandomizeTissueCoefficients(long rngSeed, bool useRngSeed, double percentage) { std::mt19937 rng; std::random_device randomDevice; if (useRngSeed) { rng.seed(rngSeed); } else { if (randomDevice.entropy() > 0.1) { rng.seed(randomDevice()); } else { rng.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count()); } } std::normal_distribution<> percentageDistribution(1, percentage / 100); for (int y = 0; y < m_TissueParameters->GetYDim(); y++) { for (int x = 0; x < m_TissueParameters->GetXDim(); x++) { for (int z = 0; z < m_TissueParameters->GetZDim(); z++) { m_AbsorptionVolume->SetData(m_AbsorptionVolume->GetData(x, y, z) * percentageDistribution(rng), x, y, z); m_ScatteringVolume->SetData(m_ScatteringVolume->GetData(x, y, z) * percentageDistribution(rng), x, y, z); } } } } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetScatteringVolume() { return m_ScatteringVolume; } mitk::pa::Volume::Pointer mitk::pa::InSilicoTissueVolume::GetAnisotropyVolume() { return m_AnisotropyVolume; } void mitk::pa::InSilicoTissueVolume::SetAbsorptionVolume(Volume::Pointer volume) { m_AbsorptionVolume = volume; } void mitk::pa::InSilicoTissueVolume::SetScatteringVolume(Volume::Pointer volume) { m_ScatteringVolume = volume; } void mitk::pa::InSilicoTissueVolume::SetAnisotropyVolume(Volume::Pointer volume) { m_AnisotropyVolume = volume; } void mitk::pa::InSilicoTissueVolume::SetSegmentationVolume(Volume::Pointer volume) { m_SegmentationVolume = volume; } diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp index 580cc4fb01..42d72b557b 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp @@ -1,182 +1,188 @@ /*=================================================================== 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" 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); } - auto generatedVolume = mitk::pa::InSilicoTissueVolume::New(parameters); - - const double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2; - 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()); } } - 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()); + auto generatedVolume = mitk::pa::InSilicoTissueVolume::New(parameters, &randomNumberGenerator); + + const double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2; + + 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_int_distribution borderTypeDistribution(0, 3); int numberOfBloodVessels = randomNumVesselDistribution(randomNumberGenerator); generatedVolume->AddIntProperty("numberOfVessels", numberOfBloodVessels); generatedVolume->AddIntProperty("bifurcationFrequency", parameters->GetVesselBifurcationFrequency()); 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); 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); /*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, 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/mitkPATissueGeneratorParameters.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPATissueGeneratorParameters.cpp index b7f8a22f78..af4359d8c2 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPATissueGeneratorParameters.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPATissueGeneratorParameters.cpp @@ -1,78 +1,79 @@ /*=================================================================== 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_BackgroundAbsorption = 0.1; + 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; } mitk::pa::TissueGeneratorParameters::~TissueGeneratorParameters() { } diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp index 5355b82627..6ca3674123 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticComposedVolumeTest.cpp @@ -1,148 +1,149 @@ /*=================================================================== 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 "mitkPAComposedVolume.h" #include "mitkIOUtil.h" #include "mitkImageReadAccessor.h" #include class mitkPhotoacousticComposedVolumeTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticComposedVolumeTestSuite); MITK_TEST(TestCreateAndDestructComposedVolume); MITK_TEST(TestAccessInvalidFluenceComponent); MITK_TEST(TestAccessInvalidFluenceComponentIndex); MITK_TEST(TestAddMultiplePairs); MITK_TEST(TestSortFunctionality); MITK_TEST(TestAccessInvalidFluenceComponentForYOffset); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::ComposedVolume::Pointer m_ComposedVolume; mitk::pa::TissueGeneratorParameters::Pointer m_DefaultParameters; mitk::pa::InSilicoTissueVolume::Pointer m_InSilicoTissueVolume; public: void setUp() { m_DefaultParameters = mitk::pa::TissueGeneratorParameters::New(); m_DefaultParameters->SetXDim(5); m_DefaultParameters->SetYDim(5); m_DefaultParameters->SetZDim(5); - m_InSilicoTissueVolume = mitk::pa::InSilicoTissueVolume::New(m_DefaultParameters); + auto rng = std::mt19937(); + m_InSilicoTissueVolume = mitk::pa::InSilicoTissueVolume::New(m_DefaultParameters, &rng); m_ComposedVolume = mitk::pa::ComposedVolume::New(m_InSilicoTissueVolume); } mitk::pa::FluenceYOffsetPair::Pointer createFluenceYOffsetPair(double value, double yOffset) { double* data = new double[125]; for (int i = 0; i < 125; ++i) data[i] = value; mitk::pa::Volume::Pointer volume = mitk::pa::Volume::New(data, 5, 5, 5, 1); return mitk::pa::FluenceYOffsetPair::New(volume, yOffset); } void TestCreateAndDestructComposedVolume() { CPPUNIT_ASSERT(m_ComposedVolume->GetNumberOfFluenceComponents() == 0); } void TestAccessInvalidFluenceComponent() { bool caughtException = false; try { m_ComposedVolume->GetFluenceValue(0, 0, 0, 0); } catch (mitk::Exception e) { caughtException = true; } CPPUNIT_ASSERT(caughtException); } void TestAddMultiplePairs() { m_ComposedVolume->AddSlice(createFluenceYOffsetPair(0, 0)); CPPUNIT_ASSERT(m_ComposedVolume->GetNumberOfFluenceComponents() == 1); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(1, 1)); CPPUNIT_ASSERT(m_ComposedVolume->GetNumberOfFluenceComponents() == 2); } void TestSortFunctionality() { m_ComposedVolume->AddSlice(createFluenceYOffsetPair(2, 2)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(-1, -1)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(1, 1)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(0, 0)); m_ComposedVolume->AddSlice(createFluenceYOffsetPair(-2, -2)); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(0, 0, 2, 0) == 2); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(1, 0, 2, 0) == -1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(2, 0, 2, 0) == 1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(3, 0, 2, 0) == 0); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(4, 0, 2, 0) == -2); m_ComposedVolume->Sort(); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(0, 0, 2, 0) == -2); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(1, 0, 2, 0) == -1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(2, 0, 2, 0) == 0); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(3, 0, 2, 0) == 1); CPPUNIT_ASSERT(m_ComposedVolume->GetFluenceValue(4, 0, 2, 0) == 2); } void TestAccessInvalidFluenceComponentIndex() { #ifdef _DEBUG m_ComposedVolume->AddSlice(createFluenceYOffsetPair(0, 0)); bool caughtException = false; try { double unusedValue = m_ComposedVolume->GetFluenceValue(0, 1, 2, 300); unusedValue = 0; } catch (mitk::Exception e) { caughtException = true; } CPPUNIT_ASSERT(caughtException); #endif } void TestAccessInvalidFluenceComponentForYOffset() { bool caughtException = false; try { m_ComposedVolume->GetYOffsetForFluenceComponentInPixels(0); } catch (mitk::Exception e) { caughtException = true; } CPPUNIT_ASSERT(caughtException); } void tearDown() { m_ComposedVolume = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticComposedVolume) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp index 5010b20d7d..0e41a51705 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp @@ -1,187 +1,189 @@ /*=================================================================== 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 class mitkPhotoacousticIOTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticIOTestSuite); MITK_TEST(testLoadInSilicoTissueNrrdFile); MITK_TEST(testLoad3DVolumeNrrdFile); MITK_TEST(testLoad3DVolumeNrrdFileWithBlur); MITK_TEST(testGetNumberOfNrrdFilesInTestDir); MITK_TEST(testGetChildFoldersFromFolder); MITK_TEST(testLoadFCMs); CPPUNIT_TEST_SUITE_END(); private: const std::string TEST_FOLDER_PATH = "testFiles/"; const std::string TEST_IN_SILICO_VOLUME_PATH = "testInSilicoVolume"; const std::string TEST_3D_Volume_PATH = "test3DVolume"; const std::string TEST_FILE_ENDING = ".nrrd"; const std::string TEST_QUALIFIED_FOLDER_PATH = TEST_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + "/"; const std::string FOLDER_FOLDER = "folder/"; const std::string FCM_PATH = TEST_FOLDER_PATH + "fcms/"; const int NUMBER_OF_NRRD_FILES_IN_TEST_DIR = 2; mitk::pa::TissueGeneratorParameters::Pointer m_VolumeProperties; mitk::pa::InSilicoTissueVolume::Pointer m_TestInSilicoVolume; mitk::pa::Volume::Pointer m_Test3DVolume; public: void setUp() { m_VolumeProperties = createTestVolumeParameters(); - m_TestInSilicoVolume = mitk::pa::InSilicoTissueVolume::New(m_VolumeProperties); + auto rng = std::mt19937(); + m_TestInSilicoVolume = mitk::pa::InSilicoTissueVolume::New(m_VolumeProperties, &rng); m_Test3DVolume = createTest3DVolume(5); itk::FileTools::CreateDirectory(TEST_FOLDER_PATH); itk::FileTools::CreateDirectory(TEST_QUALIFIED_FOLDER_PATH); itk::FileTools::CreateDirectory(TEST_FOLDER_PATH + FOLDER_FOLDER + FOLDER_FOLDER); itk::FileTools::CreateDirectory(FCM_PATH); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH)); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_QUALIFIED_FOLDER_PATH)); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH + FOLDER_FOLDER + FOLDER_FOLDER)); CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(FCM_PATH)); mitk::IOUtil::Save(m_TestInSilicoVolume->ConvertToMitkImage(), TEST_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + TEST_FILE_ENDING); mitk::IOUtil::Save(m_Test3DVolume->AsMitkImage(), TEST_FOLDER_PATH + TEST_3D_Volume_PATH + TEST_FILE_ENDING); auto yo0 = createTest3DVolume(1)->AsMitkImage(); auto yo1 = createTest3DVolume(2)->AsMitkImage(); yo0->GetPropertyList()->SetStringProperty("y-offset", "0"); yo1->GetPropertyList()->SetStringProperty("y-offset", "1"); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("y-offset")); mitk::IOUtil::Save(yo0, TEST_QUALIFIED_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + "_yo0" + TEST_FILE_ENDING); mitk::IOUtil::Save(yo1, TEST_QUALIFIED_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + "_yo1" + TEST_FILE_ENDING); } mitk::pa::Volume::Pointer createTest3DVolume(double value) { unsigned int xDim = 10; unsigned int yDim = 10; unsigned int zDim = 10; unsigned int length = xDim * yDim * zDim; double* data = new double[length]; for (unsigned int i = 0; i < length; i++) data[i] = value; return mitk::pa::Volume::New(data, xDim, yDim, zDim, 1); } mitk::pa::TissueGeneratorParameters::Pointer createTestVolumeParameters() { auto returnParameters = mitk::pa::TissueGeneratorParameters::New(); returnParameters->SetXDim(10); returnParameters->SetYDim(10); returnParameters->SetZDim(10); - returnParameters->SetBackgroundAbsorption(0); + returnParameters->SetMinBackgroundAbsorption(0); + returnParameters->SetMaxBackgroundAbsorption(0); returnParameters->SetBackgroundScattering(0); returnParameters->SetBackgroundAnisotropy(0); return returnParameters; } void assertEqual(mitk::pa::Volume::Pointer first, mitk::pa::Volume::Pointer second) { CPPUNIT_ASSERT(first->GetXDim() == second->GetXDim()); CPPUNIT_ASSERT(first->GetYDim() == second->GetYDim()); CPPUNIT_ASSERT(first->GetZDim() == second->GetZDim()); for (unsigned int x = 0; x < first->GetXDim(); ++x) for (unsigned int y = 0; y < first->GetYDim(); ++y) for (unsigned int z = 0; z < first->GetZDim(); ++z) { std::string message = "Expected " + std::to_string(first->GetData(x, y, z)) + " but was " + std::to_string(second->GetData(x, y, z)); CPPUNIT_ASSERT_MESSAGE(message, abs(first->GetData(x, y, z) - second->GetData(x, y, z)) < 1e-6); } } void testLoadInSilicoTissueNrrdFile() { auto loadedVolume = mitk::pa::IOUtil::LoadInSilicoTissueVolumeFromNrrdFile(TEST_FOLDER_PATH + TEST_IN_SILICO_VOLUME_PATH + TEST_FILE_ENDING); CPPUNIT_ASSERT(loadedVolume->GetTDim() == m_TestInSilicoVolume->GetTDim()); assertEqual(m_TestInSilicoVolume->GetAbsorptionVolume(), loadedVolume->GetAbsorptionVolume()); assertEqual(m_TestInSilicoVolume->GetScatteringVolume(), loadedVolume->GetScatteringVolume()); assertEqual(m_TestInSilicoVolume->GetAnisotropyVolume(), loadedVolume->GetAnisotropyVolume()); } void testLoad3DVolumeNrrdFile() { auto loadedVolume = mitk::pa::IOUtil::LoadNrrd(TEST_FOLDER_PATH + TEST_3D_Volume_PATH + TEST_FILE_ENDING); assertEqual(loadedVolume, m_Test3DVolume); } void testLoad3DVolumeNrrdFileWithBlur() { auto loadedVolume = mitk::pa::IOUtil::LoadNrrd(TEST_FOLDER_PATH + TEST_3D_Volume_PATH + TEST_FILE_ENDING, 1); assertEqual(loadedVolume, m_Test3DVolume); } void testGetNumberOfNrrdFilesInTestDir() { int numberOfFiles = mitk::pa::IOUtil::GetNumberOfNrrdFilesInDirectory(TEST_FOLDER_PATH); CPPUNIT_ASSERT(numberOfFiles == NUMBER_OF_NRRD_FILES_IN_TEST_DIR); } void testGetChildFoldersFromFolder() { std::vector childFolders = mitk::pa::IOUtil::GetAllChildfoldersFromFolder(TEST_FOLDER_PATH); CPPUNIT_ASSERT(childFolders.size() == 1); CPPUNIT_ASSERT(childFolders[0] == TEST_FOLDER_PATH); childFolders = mitk::pa::IOUtil::GetAllChildfoldersFromFolder(TEST_FOLDER_PATH + FOLDER_FOLDER); MITK_INFO << "ChildFolders: " << childFolders.size(); CPPUNIT_ASSERT(childFolders.size() == 1); CPPUNIT_ASSERT(childFolders[0] == TEST_FOLDER_PATH + FOLDER_FOLDER + "/folder"); } void testLoadFCMs() { auto fcm1 = createTest3DVolume(1); auto fcm2 = createTest3DVolume(2); auto fcm3 = createTest3DVolume(3); auto fcm4 = createTest3DVolume(4); mitk::IOUtil::Save(fcm1->AsMitkImage(), FCM_PATH + "fcm1_p0,0,0FluenceContributionMap.nrrd"); mitk::IOUtil::Save(fcm2->AsMitkImage(), FCM_PATH + "fcm1_p0,0,1FluenceContributionMap.nrrd"); mitk::IOUtil::Save(fcm3->AsMitkImage(), FCM_PATH + "fcm1_p1,0,0FluenceContributionMap.nrrd"); mitk::IOUtil::Save(fcm4->AsMitkImage(), FCM_PATH + "fcm1_p1,0,1FluenceContributionMap.nrrd"); int prog = 0; auto map = mitk::pa::IOUtil::LoadFluenceContributionMaps(FCM_PATH, 0, &prog); assertEqual(fcm1, map[mitk::pa::IOUtil::Position{ 0,0 }]); assertEqual(fcm2, map[mitk::pa::IOUtil::Position{ 0,1 }]); assertEqual(fcm3, map[mitk::pa::IOUtil::Position{ 1,0 }]); assertEqual(fcm4, map[mitk::pa::IOUtil::Position{ 1,1 }]); } void tearDown() { //CPPUNIT_ASSERT_MESSAGE("Resource leak of test files onto hard drive..", itksys::SystemTools::RemoveADirectory(TEST_FOLDER_PATH) == true); } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticIO) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticTissueGeneratorTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticTissueGeneratorTest.cpp index 6a1f07d09f..9606ada406 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticTissueGeneratorTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticTissueGeneratorTest.cpp @@ -1,83 +1,85 @@ /*=================================================================== 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 "mitkPATissueGenerator.h" class mitkPhotoacousticTissueGeneratorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticTissueGeneratorTestSuite); MITK_TEST(testCallWithEmptyParameters); MITK_TEST(testCallWithWorkingParameters); CPPUNIT_TEST_SUITE_END(); private: public: void setUp() { } mitk::pa::TissueGeneratorParameters::Pointer createRandomTestVolumeParameters() { auto returnParameters = mitk::pa::TissueGeneratorParameters::New(); returnParameters->SetXDim(rand() % 50 + 1); returnParameters->SetYDim(rand() % 50 + 1); returnParameters->SetZDim(rand() % 50 + 1); - returnParameters->SetBackgroundAbsorption(rand() % 100 / 10.0); + double absorb = rand() % 100 / 10.0; + returnParameters->SetMinBackgroundAbsorption(absorb); + returnParameters->SetMaxBackgroundAbsorption(absorb); returnParameters->SetBackgroundScattering(rand() % 100 / 10.0); returnParameters->SetBackgroundAnisotropy(rand() % 100 / 10.0); int min = rand() % 10; returnParameters->SetMinNumberOfVessels(min); returnParameters->SetMaxNumberOfVessels(min + (rand() % 10)); returnParameters->SetCalculateNewVesselPositionCallback( &mitk::pa::VesselMeanderStrategy::CalculateRandomlyDivergingPosition); returnParameters->SetMinVesselZOrigin(rand() % 3 + 1); returnParameters->SetMaxVesselZOrigin(rand() % 3 + 1); int minRad = rand() % 100; returnParameters->SetMinVesselRadiusInMillimeters(minRad); returnParameters->SetMaxVesselRadiusInMillimeters(minRad + (rand() % 100)); returnParameters->SetVoxelSpacingInCentimeters(1); return returnParameters; } void testCallWithEmptyParameters() { auto parameters = mitk::pa::TissueGeneratorParameters::New(); auto volume = mitk::pa::InSilicoTissueGenerator::GenerateInSilicoData(parameters); CPPUNIT_ASSERT(volume.IsNotNull()); } void testCallWithWorkingParameters() { for (int i = 0; i < 20; i++) { auto parameters = createRandomTestVolumeParameters(); auto volume = mitk::pa::InSilicoTissueGenerator::GenerateInSilicoData(parameters); CPPUNIT_ASSERT(volume.IsNotNull()); } } void tearDown() { } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticTissueGenerator) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp index 1f95d64c9b..51b0b3830c 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTest.cpp @@ -1,240 +1,242 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include "mitkPAInSilicoTissueVolume.h" #include "mitkPAVector.h" #include "mitkPAVessel.h" #include "mitkIOUtil.h" class mitkPhotoacousticVesselTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticVesselTestSuite); MITK_TEST(testEmptyInitializationProperties); MITK_TEST(testWalkInStraightLine); MITK_TEST(testBifurcate); MITK_TEST(testPartialVolume); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::Vessel::Pointer m_TestVessel; mitk::pa::Vessel::CalculateNewVesselPositionCallback m_StraightLine; mitk::pa::Vessel::CalculateNewVesselPositionCallback m_Diverging; mitk::pa::InSilicoTissueVolume::Pointer m_TestInSilicoVolume; mitk::pa::TissueGeneratorParameters::Pointer m_TestVolumeParameters; public: void setUp() { auto params = mitk::pa::VesselProperties::New(); m_TestVessel = mitk::pa::Vessel::New(params); m_StraightLine = &mitk::pa::VesselMeanderStrategy::CalculateNewPositionInStraightLine; m_Diverging = &mitk::pa::VesselMeanderStrategy::CalculateRandomlyDivergingPosition; m_TestVolumeParameters = createTestVolumeParameters(); - m_TestInSilicoVolume = mitk::pa::InSilicoTissueVolume::New(m_TestVolumeParameters); + auto rng = std::mt19937(); + m_TestInSilicoVolume = mitk::pa::InSilicoTissueVolume::New(m_TestVolumeParameters, &rng); } mitk::pa::TissueGeneratorParameters::Pointer createTestVolumeParameters() { auto returnParameters = mitk::pa::TissueGeneratorParameters::New(); returnParameters->SetXDim(10); returnParameters->SetYDim(10); returnParameters->SetZDim(10); - returnParameters->SetBackgroundAbsorption(0); + returnParameters->SetMinBackgroundAbsorption(0); + returnParameters->SetMaxBackgroundAbsorption(0); returnParameters->SetBackgroundScattering(0); returnParameters->SetBackgroundAnisotropy(0); return returnParameters; } void testEmptyInitializationProperties() { CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == true); } void testPartialVolume() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); params->SetDoPartialVolume(true); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(100); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 3)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 3) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 5) - 5.85786438) <= 1e-6); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 3)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 3) - 5.85786438) <= 1e-6); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 3)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 3) - 5.85786438) <= 1e-6); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 3, 5) - 5.85786438) <= 1e-6); } void testWalkInStraightLine() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); params->SetDoPartialVolume(false); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(100); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); mitk::IOUtil::Save(m_TestInSilicoVolume->GetAbsorptionVolume()->AsMitkImage(), "C:/Users/groehl/Desktop/test.nrrd"); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 5, 4)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 6, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 6, 4)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 5)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 6)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 6)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 5, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 5, 4)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 6, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 6, 4)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 5)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 5)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 6)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 6)) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(1, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(2, 4, 4) - 10) <= mitk::eps); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(3, 4, 4)), abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(3, 4, 4)) <= mitk::eps); } void testBifurcate() { auto testPosition = mitk::pa::Vector::New(); testPosition->SetElement(0, 0); testPosition->SetElement(1, 4); testPosition->SetElement(2, 4); auto testDirection = mitk::pa::Vector::New(); testDirection->SetElement(0, 1); testDirection->SetElement(1, 0); testDirection->SetElement(2, 0); auto params = mitk::pa::VesselProperties::New(); params->SetRadiusInVoxel(1); params->SetBifurcationFrequency(1); params->SetAbsorptionCoefficient(10); params->SetScatteringCoefficient(10); params->SetAnisotopyCoefficient(10); params->SetPositionVector(testPosition); params->SetDirectionVector(testDirection); m_TestVessel = mitk::pa::Vessel::New(params); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(m_TestVessel->IsFinished() == false); CPPUNIT_ASSERT(abs(m_TestInSilicoVolume->GetAbsorptionVolume()->GetData(0, 4, 4)) <= mitk::eps); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); m_TestVessel->ExpandVessel(m_TestInSilicoVolume, m_StraightLine, 0, nullptr); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == true); std::mt19937 rng; auto bifurcationVessel = m_TestVessel->Bifurcate(&rng); CPPUNIT_ASSERT(m_TestVessel->CanBifurcate() == false); CPPUNIT_ASSERT(bifurcationVessel->CanBifurcate() == false); } void tearDown() { } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVessel) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTreeTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTreeTest.cpp index 96884e36b7..de8cde4af0 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTreeTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVesselTreeTest.cpp @@ -1,106 +1,108 @@ /*=================================================================== 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 "mitkPAVesselTree.h" #include "mitkPAInSilicoTissueVolume.h" using namespace mitk::pa; class mitkPhotoacousticVesselTreeTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticVesselTreeTestSuite); MITK_TEST(testVesselTreeInitialBehavior); MITK_TEST(testCallStepMethod); CPPUNIT_TEST_SUITE_END(); private: public: VesselTree::Pointer m_Tree; VesselProperties::Pointer m_VesselProperties; Vessel::CalculateNewVesselPositionCallback m_StraightLine; InSilicoTissueVolume::Pointer m_TestInSilicoVolume; void setUp() { m_VesselProperties = VesselProperties::New(); m_Tree = VesselTree::New(m_VesselProperties); m_StraightLine = &VesselMeanderStrategy::CalculateNewPositionInStraightLine; - m_TestInSilicoVolume = InSilicoTissueVolume::New(createTestVolumeParameters()); + auto rng = std::mt19937(); + m_TestInSilicoVolume = InSilicoTissueVolume::New(createTestVolumeParameters(), &rng); } TissueGeneratorParameters::Pointer createTestVolumeParameters() { auto returnParameters = TissueGeneratorParameters::New(); returnParameters->SetXDim(10); returnParameters->SetYDim(10); returnParameters->SetZDim(10); - returnParameters->SetBackgroundAbsorption(0); + returnParameters->SetMinBackgroundAbsorption(0); + returnParameters->SetMaxBackgroundAbsorption(0); returnParameters->SetBackgroundScattering(0); returnParameters->SetBackgroundAnisotropy(0); return returnParameters; } void testVesselTreeInitialBehavior() { CPPUNIT_ASSERT(m_Tree->IsFinished() == true); } void testCallStepMethod() { //really bad test atm.. The only thing that is tested is that the method can be called without a crash. //But hey - it is something :P m_VesselProperties->SetRadiusInVoxel(2); std::mt19937 rng; rng.seed(1); m_Tree = VesselTree::New(m_VesselProperties); m_Tree->Step(m_TestInSilicoVolume, m_StraightLine, 0, &rng); rng.seed(1); auto secondTree = VesselTree::New(m_VesselProperties); secondTree->Step(m_TestInSilicoVolume, m_StraightLine, 0, &rng); CPPUNIT_ASSERT(Equal(m_Tree, secondTree, 1e-6, true)); secondTree->Step(m_TestInSilicoVolume, m_StraightLine, 0, &rng); CPPUNIT_ASSERT(!Equal(m_Tree, secondTree, 1e-6, true)); int i = 0; for (; i < 1000; i++) { secondTree->Step(m_TestInSilicoVolume, m_StraightLine, 0, &rng); if (secondTree->IsFinished()) break; } CPPUNIT_ASSERT(i < 999); } void tearDown() { } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVesselTree) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp index 36e086dff3..d5dc884b51 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVolumeTest.cpp @@ -1,347 +1,385 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include "mitkPAInSilicoTissueVolume.h" #include "mitkPATissueGeneratorParameters.h" class mitkPhotoacousticVolumeTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticVolumeTestSuite); + MITK_TEST(TestUniformDistributionIsUniform); MITK_TEST(TestInitializedTissueContainsOnlyZeros); MITK_TEST(TestConvertedMitkImageContainsOnlyZerosOrAir); MITK_TEST(TestTissueVolumeContainsCorrectAbsorptionNumber); MITK_TEST(TestTissueVolumeContainsCorrectScatteringNumber); MITK_TEST(TestTissueVolumeContainsCorrectAnisotropyNumber); MITK_TEST(testSecondConstructor); MITK_TEST(testCompleteAirVoxelInclusion); MITK_TEST(testHalfAirVoxelInclusion); MITK_TEST(testCompleteAirAndSkinVoxelInclusion); MITK_TEST(testRandomizeCoefficients); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::InSilicoTissueVolume::Pointer m_PhotoacousticVolume; mitk::pa::TissueGeneratorParameters::Pointer m_TissueGeneratorParameters; public: void setUp() { m_TissueGeneratorParameters = mitk::pa::TissueGeneratorParameters::New(); - m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); + auto rng = std::mt19937(); + m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters, &rng); + } + + void TestUniformDistributionIsUniform() + { + + int dims = 30; + m_TissueGeneratorParameters->SetXDim(dims); + m_TissueGeneratorParameters->SetYDim(dims); + m_TissueGeneratorParameters->SetZDim(dims); + m_TissueGeneratorParameters->SetAirThicknessInMillimeters(0); + + m_TissueGeneratorParameters->SetMinBackgroundAbsorption(0.001); + m_TissueGeneratorParameters->SetMaxBackgroundAbsorption(0.2); + + auto rng = std::mt19937(); + m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters, &rng); + + for (int x = 0; x < dims; x++) + { + for (int y = 0; y < dims; y++) + { + for (int z = 0; z < dims; z++) + { + CPPUNIT_ASSERT_EQUAL_MESSAGE("Every absorption should be in bounds.", + m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(x, y, z) >= 0.001 && + m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(x, y, z) <= 0.2, true); + } + } + } + } void TestInitializedTissueContainsOnlyZeros() { int dims = 30; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(0); - m_TissueGeneratorParameters->SetBackgroundAbsorption(0); - m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); + m_TissueGeneratorParameters->SetMinBackgroundAbsorption(0); + m_TissueGeneratorParameters->SetMaxBackgroundAbsorption(0); + auto rng = std::mt19937(); + m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters, &rng); for (int x = 0; x < dims; x++) { for (int y = 0; y < dims; y++) { for (int z = 0; z < dims; z++) { CPPUNIT_ASSERT_EQUAL_MESSAGE("Every field should be initialized with 0.", abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(x, y, z)) < mitk::eps, true); } } } } void TestConvertedMitkImageContainsOnlyZerosOrAir() { int dims = 30; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); - m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); + auto rng = std::mt19937(); + m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters, &rng); mitk::Image::Pointer testImage = m_PhotoacousticVolume->ConvertToMitkImage(); mitk::ImageReadAccessor imgMemAcc(testImage); double* imagePointer = (double*)imgMemAcc.GetData(); for (int index = 0; index < dims*dims*dims; index++, imagePointer++) { CPPUNIT_ASSERT_EQUAL_MESSAGE("Every voxel in image should be 0.1 or 0.0001.", true, abs(*imagePointer - 0.1) <= mitk::eps || abs(*imagePointer - 0.0001) <= mitk::eps); } } void TestTissueVolumeContainsCorrectAbsorptionNumber() { int dims = 2; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); - m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); + auto rng = std::mt19937(); + m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters, &rng); m_PhotoacousticVolume->SetVolumeValues(0, 0, 0, 0, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 0, 1, 1, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 0, 2, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 1, 3, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 0, 4, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 1, 5, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 0, 6, 0, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 1, 7, 0, 0); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 0.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 1.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 2.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 3.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 1, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 4.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 5.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 6.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 7.0, m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1)); } void TestTissueVolumeContainsCorrectScatteringNumber() { int dims = 2; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); - m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); + auto rng = std::mt19937(); + m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters, &rng); m_PhotoacousticVolume->SetVolumeValues(0, 0, 0, 0, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 0, 1, 0, 1, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 0, 0, 2, 0); m_PhotoacousticVolume->SetVolumeValues(0, 1, 1, 0, 3, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 0, 0, 4, 0); m_PhotoacousticVolume->SetVolumeValues(1, 0, 1, 0, 5, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 0, 0, 6, 0); m_PhotoacousticVolume->SetVolumeValues(1, 1, 1, 0, 7, 0); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 0.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 1.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 2.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 3.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 1, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 4.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 5.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 6.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 7.0, m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1)); } void TestTissueVolumeContainsCorrectAnisotropyNumber() { int dims = 2; m_TissueGeneratorParameters->SetXDim(dims); m_TissueGeneratorParameters->SetYDim(dims); m_TissueGeneratorParameters->SetZDim(dims); - m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters); + auto rng = std::mt19937(); + m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(m_TissueGeneratorParameters, &rng); m_PhotoacousticVolume->SetVolumeValues(0, 0, 0, 0, 0, 0); m_PhotoacousticVolume->SetVolumeValues(0, 0, 1, 0, 0, 1); m_PhotoacousticVolume->SetVolumeValues(0, 1, 0, 0, 0, 2); m_PhotoacousticVolume->SetVolumeValues(0, 1, 1, 0, 0, 3); m_PhotoacousticVolume->SetVolumeValues(1, 0, 0, 0, 0, 4); m_PhotoacousticVolume->SetVolumeValues(1, 0, 1, 0, 0, 5); m_PhotoacousticVolume->SetVolumeValues(1, 1, 0, 0, 0, 6); m_PhotoacousticVolume->SetVolumeValues(1, 1, 1, 0, 0, 7); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 0.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 1.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 2.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 3.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 1, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 4.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 0, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 5.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 0, 1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 6.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("Should be correct value", 7.0, m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1)); } mitk::pa::Volume::Pointer createTestVolume(double value) { double* data = new double[27]; for (int i = 0; i < 27; ++i) data[i] = value; return mitk::pa::Volume::New(data, 3, 3, 3, 1); } void assertEqual(mitk::pa::Volume::Pointer first, mitk::pa::Volume::Pointer second) { CPPUNIT_ASSERT(first->GetXDim() == second->GetXDim()); CPPUNIT_ASSERT(first->GetYDim() == second->GetYDim()); CPPUNIT_ASSERT(first->GetZDim() == second->GetZDim()); for (unsigned int x = 0; x < first->GetXDim(); ++x) for (unsigned int y = 0; y < first->GetYDim(); ++y) for (unsigned int z = 0; z < first->GetZDim(); ++z) CPPUNIT_ASSERT(abs(first->GetData(x, y, z) - second->GetData(x, y, z)) < mitk::eps); } void testSecondConstructor() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); assertEqual(m_PhotoacousticVolume->GetAbsorptionVolume(), absorption); assertEqual(m_PhotoacousticVolume->GetScatteringVolume(), scattering); assertEqual(m_PhotoacousticVolume->GetAnisotropyVolume(), anisotropy); assertEqual(m_PhotoacousticVolume->GetSegmentationVolume(), segmentation); } void testCompleteAirVoxelInclusion() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(10); m_TissueGeneratorParameters->SetSkinThicknessInMillimeters(0); m_TissueGeneratorParameters->SetAirAbsorption(2); m_TissueGeneratorParameters->SetAirScattering(4); m_TissueGeneratorParameters->SetAirAnisotropy(6); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 6) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 3) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 3) < mitk::eps); } void testRandomizeCoefficients() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(1); mitk::pa::Volume::Pointer anisotropy = createTestVolume(1); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetRandomizePhysicalProperties(true); m_TissueGeneratorParameters->SetRandomizePhysicalPropertiesPercentage(1); m_TissueGeneratorParameters->SetRngSeed(17); m_TissueGeneratorParameters->SetUseRngSeed(true); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 1) < 0.1); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 1) < 0.1); } void testCompleteAirAndSkinVoxelInclusion() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(10); m_TissueGeneratorParameters->SetSkinThicknessInMillimeters(10); m_TissueGeneratorParameters->SetAirAbsorption(2); m_TissueGeneratorParameters->SetAirScattering(4); m_TissueGeneratorParameters->SetAirAnisotropy(6); m_TissueGeneratorParameters->SetSkinAbsorption(4); m_TissueGeneratorParameters->SetSkinScattering(8); m_TissueGeneratorParameters->SetSkinAnisotropy(12); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 8) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 6) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 12) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 3) < mitk::eps); } void testHalfAirVoxelInclusion() { mitk::pa::Volume::Pointer absorption = createTestVolume(1); mitk::pa::Volume::Pointer scattering = createTestVolume(2); mitk::pa::Volume::Pointer anisotropy = createTestVolume(3); mitk::pa::Volume::Pointer segmentation = createTestVolume(4); mitk::PropertyList::Pointer properties = mitk::PropertyList::New(); m_TissueGeneratorParameters->SetXDim(3); m_TissueGeneratorParameters->SetYDim(3); m_TissueGeneratorParameters->SetZDim(3); m_TissueGeneratorParameters->SetAirThicknessInMillimeters(15); m_TissueGeneratorParameters->SetSkinThicknessInMillimeters(0); m_TissueGeneratorParameters->SetAirAbsorption(2); m_TissueGeneratorParameters->SetAirScattering(4); m_TissueGeneratorParameters->SetAirAnisotropy(6); m_PhotoacousticVolume = mitk::pa::InSilicoTissueVolume::New(absorption, scattering, anisotropy, segmentation, m_TissueGeneratorParameters, properties); m_PhotoacousticVolume->FinalizeVolume(); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(0, 0, 0) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 1) - 1.5) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAbsorptionVolume()->GetData(1, 1, 2) - 1) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(0, 0, 0) - 4) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 1) - 3) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetScatteringVolume()->GetData(1, 1, 2) - 2) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(0, 0, 0) - 6) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 1) - 4.5) < mitk::eps); CPPUNIT_ASSERT(abs(m_PhotoacousticVolume->GetAnisotropyVolume()->GetData(1, 1, 2) - 3) < mitk::eps); } void tearDown() { m_PhotoacousticVolume = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVolume) diff --git a/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp b/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp index 66e5341de6..f611ad9a8b 100644 --- a/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkSlicedVolumeGeneratorTest.cpp @@ -1,187 +1,188 @@ ///*=================================================================== //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 class mitkSlicedVolumeGeneratorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkSlicedVolumeGeneratorTestSuite); MITK_TEST(testConstructorDestructor); MITK_TEST(testGetSlicedFluenceVolume); MITK_TEST(testGetSlicedFluenceVolumeInverse); MITK_TEST(testGetSlicedFluenceVolumeWithPrecorrection); MITK_TEST(testGetSlicedFluenceVolumeWithPrecorrectionInverse); MITK_TEST(testGetSlicedSignalVolume); MITK_TEST(testGetSlicedAbsorptionVolume); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::ComposedVolume::Pointer m_ComposedVolume; mitk::pa::TissueGeneratorParameters::Pointer m_DefaultParameters; mitk::pa::InSilicoTissueVolume::Pointer m_InSilicoTissueVolume; mitk::pa::SlicedVolumeGenerator::Pointer m_SlicedVolumeGenerator; mitk::pa::Volume::Pointer m_PrecorrectionVolume; public: void setUp() { m_SlicedVolumeGenerator = nullptr; m_DefaultParameters = mitk::pa::TissueGeneratorParameters::New(); m_DefaultParameters->SetXDim(3); m_DefaultParameters->SetYDim(3); m_DefaultParameters->SetZDim(3); - m_InSilicoTissueVolume = mitk::pa::InSilicoTissueVolume::New(m_DefaultParameters); + auto rng = std::mt19937(); + m_InSilicoTissueVolume = mitk::pa::InSilicoTissueVolume::New(m_DefaultParameters, &rng); m_ComposedVolume = mitk::pa::ComposedVolume::New(m_InSilicoTissueVolume); m_ComposedVolume->AddSlice(CreateValidationPair(-1, 1)); m_ComposedVolume->AddSlice(CreateValidationPair(0, 3)); m_ComposedVolume->AddSlice(CreateValidationPair(1, 6)); m_PrecorrectionVolume = CreatePrecorrectionVolume(); } mitk::pa::Volume::Pointer CreatePrecorrectionVolume() { double* data = new double[27]; for (int i = 0; i < 27; ++i) data[i] = 0.5; return mitk::pa::Volume::New(data, 3, 3, 3, 1); } void FillYSliceWith(mitk::pa::Volume::Pointer fluenceVolume, double ySlice, double value) { for (unsigned int x = 0; x < fluenceVolume->GetXDim(); ++x) for (unsigned int z = 0; z < fluenceVolume->GetZDim(); ++z) { fluenceVolume->SetData(value, x, ySlice, z); } } mitk::pa::FluenceYOffsetPair::Pointer CreateValidationPair(double yOffset, int start) { double* data = new double[27]; mitk::pa::Volume::Pointer fluenceVolume = mitk::pa::Volume::New(data, 3, 3, 3, 1); FillYSliceWith(fluenceVolume, 0, start + 0); FillYSliceWith(fluenceVolume, 1, start + 1); FillYSliceWith(fluenceVolume, 2, start + 2); return mitk::pa::FluenceYOffsetPair::New(fluenceVolume, yOffset); } void AssertYSliceValue(mitk::pa::Volume::Pointer fluenceVolume, double ySlice, double value) { for (unsigned int x = 0; x < fluenceVolume->GetXDim(); ++x) for (unsigned int z = 0; z < fluenceVolume->GetZDim(); ++z) { std::string msg = "Expected: " + std::to_string(value) + " actual: " + std::to_string(fluenceVolume->GetData(x, ySlice, z)); CPPUNIT_ASSERT_MESSAGE(msg, abs(fluenceVolume->GetData(x, ySlice, z) - value) < mitk::eps); } } void testConstructorDestructor() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(0, false, nullptr, false); CPPUNIT_ASSERT(m_SlicedVolumeGenerator.IsNotNull()); } void testGetSlicedFluenceVolume() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 1); AssertYSliceValue(slicedFluence, 1, 4); AssertYSliceValue(slicedFluence, 2, 8); } void testGetSlicedFluenceVolumeInverse() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, true); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 1); AssertYSliceValue(slicedFluence, 1, 1.0 / 4.0); AssertYSliceValue(slicedFluence, 2, 1.0 / 8.0); } void testGetSlicedFluenceVolumeWithPrecorrection() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, true, m_PrecorrectionVolume, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 2); AssertYSliceValue(slicedFluence, 1, 8); AssertYSliceValue(slicedFluence, 2, 16); } void testGetSlicedFluenceVolumeWithPrecorrectionInverse() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, true, m_PrecorrectionVolume, true); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedFluenceImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); AssertYSliceValue(slicedFluence, 0, 1.0 / 2); AssertYSliceValue(slicedFluence, 1, 1.0 / 8); AssertYSliceValue(slicedFluence, 2, 1.0 / 16); } void testGetSlicedSignalVolume() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedSignalImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); - AssertYSliceValue(slicedFluence, 0, 1 * m_DefaultParameters->GetBackgroundAbsorption()); - AssertYSliceValue(slicedFluence, 1, 4 * m_DefaultParameters->GetBackgroundAbsorption()); - AssertYSliceValue(slicedFluence, 2, 8 * m_DefaultParameters->GetBackgroundAbsorption()); + AssertYSliceValue(slicedFluence, 0, 1 * m_DefaultParameters->GetMinBackgroundAbsorption()); + AssertYSliceValue(slicedFluence, 1, 4 * m_DefaultParameters->GetMinBackgroundAbsorption()); + AssertYSliceValue(slicedFluence, 2, 8 * m_DefaultParameters->GetMinBackgroundAbsorption()); } void testGetSlicedAbsorptionVolume() { m_SlicedVolumeGenerator = mitk::pa::SlicedVolumeGenerator::New(1, false, nullptr, false); mitk::pa::Volume::Pointer slicedFluence = m_SlicedVolumeGenerator->GetSlicedGroundTruthImageFromComposedVolume(m_ComposedVolume); CPPUNIT_ASSERT(slicedFluence->GetXDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetYDim() == 3); CPPUNIT_ASSERT(slicedFluence->GetZDim() == 3); - AssertYSliceValue(slicedFluence, 0, m_DefaultParameters->GetBackgroundAbsorption()); - AssertYSliceValue(slicedFluence, 1, m_DefaultParameters->GetBackgroundAbsorption()); - AssertYSliceValue(slicedFluence, 2, m_DefaultParameters->GetBackgroundAbsorption()); + AssertYSliceValue(slicedFluence, 0, m_DefaultParameters->GetMinBackgroundAbsorption()); + AssertYSliceValue(slicedFluence, 1, m_DefaultParameters->GetMinBackgroundAbsorption()); + AssertYSliceValue(slicedFluence, 2, m_DefaultParameters->GetMinBackgroundAbsorption()); } void tearDown() { m_SlicedVolumeGenerator = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkSlicedVolumeGenerator)