diff --git a/Modules/PhotoacousticsLib/CMakeLists.txt b/Modules/PhotoacousticsLib/CMakeLists.txt index e40d0232c7..8c89b2658e 100644 --- a/Modules/PhotoacousticsLib/CMakeLists.txt +++ b/Modules/PhotoacousticsLib/CMakeLists.txt @@ -1,13 +1,12 @@ MITK_CREATE_MODULE( INCLUDE_DIRS PUBLIC include INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} DEPENDS PUBLIC MitkAlgorithmsExt tinyxml PACKAGE_DEPENDS tinyxml PUBLIC ITK - PUBLIC Boost ) -add_subdirectory(MCxyz) +add_subdirectory(MitkMCxyz) add_subdirectory(test) diff --git a/Modules/PhotoacousticsLib/MCxyz/CMakeLists.txt b/Modules/PhotoacousticsLib/MitkMCxyz/CMakeLists.txt similarity index 92% rename from Modules/PhotoacousticsLib/MCxyz/CMakeLists.txt rename to Modules/PhotoacousticsLib/MitkMCxyz/CMakeLists.txt index 07002f170d..fdbf87e92e 100644 --- a/Modules/PhotoacousticsLib/MCxyz/CMakeLists.txt +++ b/Modules/PhotoacousticsLib/MitkMCxyz/CMakeLists.txt @@ -1,11 +1,11 @@ OPTION(BUILD_PhotoacousticSimulationMCxyz "Build MiniApp for Photoacoustic Simulation Module (mcxyz)" ON) IF(BUILD_PhotoacousticSimulationMCxyz) PROJECT( MitkMCxyz ) mitk_create_executable(MCxyz DEPENDS MitkCommandLine MitkCore MitkPhotoacousticsLib PACKAGE_DEPENDS - CPP_FILES MCxyz.cpp) + CPP_FILES MitkMCxyz.cpp) install(TARGETS ${EXECUTABLE_TARGET} RUNTIME DESTINATION bin) ENDIF() diff --git a/Modules/PhotoacousticsLib/MCxyz/MCxyz.cpp b/Modules/PhotoacousticsLib/MitkMCxyz/MitkMCxyz.cpp similarity index 99% rename from Modules/PhotoacousticsLib/MCxyz/MCxyz.cpp rename to Modules/PhotoacousticsLib/MitkMCxyz/MitkMCxyz.cpp index a40f1ee24a..fc2fe8b125 100644 --- a/Modules/PhotoacousticsLib/MCxyz/MCxyz.cpp +++ b/Modules/PhotoacousticsLib/MitkMCxyz/MitkMCxyz.cpp @@ -1,1472 +1,1472 @@ /*=================================================================== 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. ===================================================================*/ // Please retain the following copyright notice /****************************************************************** * based on mcxyz.c Oct2014 * * mcxyz.c, in ANSI Standard C programing language * * created 2010, 2012 by * Steven L. JACQUES * Ting LI * Oregon Health & Science University * *******************************************************************/ #include #include #include #include #include #include #include #include #include #include #include "mitkCommandLineParser.h" #include "mitkIOUtil.h" #include "mitkImageCast.h" #include #include #include #include #include #include #include #include #ifdef __linux__ #include #include #else #include #endif #define ls 1.0E-7 /* Moving photon a little bit off the voxel face */ #define PI 3.1415926 #define ALIVE 1 /* if photon not yet terminated */ #define DEAD 0 /* if photon is to be terminated */ #define THRESHOLD 0.01 /* used in roulette */ #define CHANCE 0.1 /* used in roulette */ #define SQR(x) (x*x) #define SIGN(x) ((x)>=0 ? 1:-1) #define ONE_MINUS_COSZERO 1.0E-12 /* If 1-cos(theta) <= ONE_MINUS_COSZERO, fabs(theta) <= 1e-6 rad. */ /* If 1+cos(theta) <= ONE_MINUS_COSZERO, fabs(PI-theta) <= 1e-6 rad. */ /* Struct for storing x,y and z coordinates */ struct Location { int x, y, z; double absorb; }; struct Location initLocation(int x, int y, int z, double absorb) { struct Location loc; loc.x = x; loc.y = y; loc.z = z; loc.absorb = absorb; return loc; } class DetectorVoxel { public: Location location; std::vector* recordedPhotonRoute = new std::vector(); double* fluenceContribution; double m_PhotonNormalizationValue; long m_NumberPhotonsCurrent; DetectorVoxel(Location location, long totalNumberOfVoxels, double photonNormalizationValue) { this->location = location; this->fluenceContribution = (double *)malloc(totalNumberOfVoxels * sizeof(double)); for (int j = 0; j < totalNumberOfVoxels; j++) fluenceContribution[j] = 0; // ensure fluenceContribution[] starts empty. m_NumberPhotonsCurrent = 0; m_PhotonNormalizationValue = photonNormalizationValue; } }; bool verbose(false); class InputValues { private: std::string inputFilename; int tissueIterator; long long ix, iy, iz; public: int mcflag, launchflag, boundaryflag; double xfocus, yfocus, zfocus; double ux0, uy0, uz0; double radius; double waist; double xs, ys, zs; /* launch position */ int Nx, Ny, Nz, numberOfTissueTypes; /* # of bins */ char* tissueType; double* muaVector; double* musVector; double* gVector; double* normalizationVector; double xSpacing, ySpacing, zSpacing; double simulationTimeFromFile; long long Nphotons; long totalNumberOfVoxels; double* totalFluence; std::string myname; DetectorVoxel* detectorVoxel; mitk::Image::Pointer m_inputImage; mitk::Image::Pointer m_normalizationImage; InputValues() { tissueType = nullptr; muaVector = nullptr; musVector = nullptr; gVector = nullptr; detectorVoxel = nullptr; normalizationVector = nullptr; mcflag = 0; launchflag = 0; boundaryflag = 0; } double GetNormalizationValue(int x, int y, int z) { if (normalizationVector) return normalizationVector[z*Ny*Nx + x*Ny + y]; else return 1; } void LoadValues(std::string localInputFilename, float yOffset, std::string normalizationFilename, bool simulatePVFC) { inputFilename = localInputFilename; try { m_inputImage = mitk::IOUtil::LoadImage(inputFilename); } catch (...) { if (verbose) std::cout << "No .nrrd file found ... switching to legacy mode." << std::endl; } try { if (simulatePVFC && !normalizationFilename.empty()) m_normalizationImage = mitk::IOUtil::LoadImage(normalizationFilename); } catch (...) { if (verbose) std::cout << "No normalization .nrrd file found ... will not normalize PVFC" << std::endl; } if (m_normalizationImage.IsNotNull()) { mitk::ImageReadAccessor readAccess3(m_normalizationImage, m_normalizationImage->GetVolumeData(0)); normalizationVector = (double *)readAccess3.GetData(); } if (m_inputImage.IsNotNull()) // load stuff from nrrd file { simulationTimeFromFile = 0; Nx = m_inputImage->GetDimensions()[1]; Ny = m_inputImage->GetDimensions()[0]; Nz = m_inputImage->GetDimensions()[2]; xSpacing = m_inputImage->GetGeometry(0)->GetSpacing()[0]; ySpacing = m_inputImage->GetGeometry(0)->GetSpacing()[1]; zSpacing = m_inputImage->GetGeometry(0)->GetSpacing()[2]; mcflag = std::stoi(m_inputImage->GetProperty("mcflag")->GetValueAsString().c_str()); // mcflag, 0 = uniform, 1 = Gaussian, 2 = iso-pt, 4 = monospectral fraunhofer setup launchflag = std::stoi(m_inputImage->GetProperty("launchflag")->GetValueAsString().c_str());// 0 = let mcxyz calculate trajectory, 1 = manually set launch vector boundaryflag = std::stoi(m_inputImage->GetProperty("boundaryflag")->GetValueAsString().c_str());// 0 = no boundaries, 1 = escape at boundaries, 2 = escape at surface only xs = std::stod(m_inputImage->GetProperty("launchPointX")->GetValueAsString().c_str()); ys = std::stod(m_inputImage->GetProperty("launchPointY")->GetValueAsString().c_str()) + yOffset; zs = std::stod(m_inputImage->GetProperty("launchPointZ")->GetValueAsString().c_str()); xfocus = std::stod(m_inputImage->GetProperty("focusPointX")->GetValueAsString().c_str()); yfocus = std::stod(m_inputImage->GetProperty("focusPointY")->GetValueAsString().c_str()); zfocus = std::stod(m_inputImage->GetProperty("focusPointZ")->GetValueAsString().c_str()); ux0 = std::stod(m_inputImage->GetProperty("trajectoryVectorX")->GetValueAsString().c_str()); uy0 = std::stod(m_inputImage->GetProperty("trajectoryVectorY")->GetValueAsString().c_str()); uz0 = std::stod(m_inputImage->GetProperty("trajectoryVectorZ")->GetValueAsString().c_str()); radius = std::stod(m_inputImage->GetProperty("radius")->GetValueAsString().c_str()); waist = std::stod(m_inputImage->GetProperty("waist")->GetValueAsString().c_str()); totalNumberOfVoxels = Nx*Ny*Nz; if (verbose) std::cout << totalNumberOfVoxels << " = sizeof totalNumberOfVoxels" << std::endl; muaVector = (double *)malloc(totalNumberOfVoxels * sizeof(double)); /* tissue structure */ musVector = (double *)malloc(totalNumberOfVoxels * sizeof(double)); /* tissue structure */ gVector = (double *)malloc(totalNumberOfVoxels * sizeof(double)); /* tissue structure */ mitk::ImageReadAccessor readAccess0(m_inputImage, m_inputImage->GetVolumeData(0)); muaVector = (double *)readAccess0.GetData(); mitk::ImageReadAccessor readAccess1(m_inputImage, m_inputImage->GetVolumeData(1)); musVector = (double *)readAccess1.GetData(); mitk::ImageReadAccessor readAccess2(m_inputImage, m_inputImage->GetVolumeData(2)); gVector = (double *)readAccess2.GetData(); } else { mitkThrow() << "No longer support loading of binary tissue files."; } } }; class ReturnValues { private: long i1 = 0, i2 = 31; // used Random Generator long ma[56]; // used Random Generator /* ma[0] is not used. */ long mj, mk; short i, ii; public: long long Nphotons; double* totalFluence; std::string myname; DetectorVoxel* detectorVoxel; ReturnValues() { detectorVoxel = nullptr; Nphotons = 0; totalFluence = nullptr; } /* SUBROUTINES */ /************************************************************************** * RandomGen * A random number generator that generates uniformly * distributed random numbers between 0 and 1 inclusive. * The algorithm is based on: * W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. * Flannery, "Numerical Recipes in C," Cambridge University * Press, 2nd edition, (1992). * and * D.E. Knuth, "Seminumerical Algorithms," 2nd edition, vol. 2 * of "The Art of Computer Programming", Addison-Wesley, (1981). * * When Type is 0, sets Seed as the seed. Make sure 0 b) m = a; else m = b; return m; } /*********************************************************** * min2 ****/ double min2(double a, double b) { double m; if (a >= b) m = b; else m = a; return m; } /*********************************************************** * min3 ****/ double min3(double a, double b, double c) { double m; if (a <= min2(b, c)) m = a; else if (b <= min2(a, c)) m = b; else m = c; return m; } /******************** * my version of FindVoxelFace for no scattering. * s = ls + FindVoxelFace2(x,y,z, tempx, tempy, tempz, dx, dy, dz, ux, uy, uz); ****/ - double FindVoxelFace2(double x1, double y1, double z1, double x2, double y2, double z2, double dx, double dy, double dz, double ux, double uy, double uz) + double FindVoxelFace2(double x1, double y1, double z1, double /*x2*/, double /*y2*/, double /*z2*/, double dx, double dy, double dz, double ux, double uy, double uz) { int ix1 = floor(x1 / dx); int iy1 = floor(y1 / dy); int iz1 = floor(z1 / dz); int ix2, iy2, iz2; if (ux >= 0) ix2 = ix1 + 1; else ix2 = ix1; if (uy >= 0) iy2 = iy1 + 1; else iy2 = iy1; if (uz >= 0) iz2 = iz1 + 1; else iz2 = iz1; double xs = fabs((ix2*dx - x1) / ux); double ys = fabs((iy2*dy - y1) / uy); double zs = fabs((iz2*dz - z1) / uz); double s = min3(xs, ys, zs); return s; } /*********************************************************** * FRESNEL REFLECTANCE * Computes reflectance as photon passes from medium 1 to * medium 2 with refractive indices n1,n2. Incident * angle a1 is specified by cosine value ca1 = cos(a1). * Program returns value of transmitted angle a1 as * value in *ca2_Ptr = cos(a2). ****/ double RFresnel(double n1, /* incident refractive index.*/ double n2, /* transmit refractive index.*/ double ca1, /* cosine of the incident */ /* angle a1, 00. */ { double r; if (n1 == n2) { /** matched boundary. **/ *ca2_Ptr = ca1; r = 0.0; } else if (ca1 > (1.0 - 1.0e-12)) { /** normal incidence. **/ *ca2_Ptr = ca1; r = (n2 - n1) / (n2 + n1); r *= r; } else if (ca1 < 1.0e-6) { /** very slanted. **/ *ca2_Ptr = 0.0; r = 1.0; } else { /** general. **/ double sa1, sa2; /* sine of incident and transmission angles. */ double ca2; /* cosine of transmission angle. */ sa1 = sqrt(1 - ca1*ca1); sa2 = n1*sa1 / n2; if (sa2 >= 1.0) { /* double check for total internal reflection. */ *ca2_Ptr = 0.0; r = 1.0; } else { double cap, cam; /* cosines of sum ap or diff am of the two */ /* angles: ap = a1 + a2, am = a1 - a2. */ double sap, sam; /* sines. */ *ca2_Ptr = ca2 = sqrt(1 - sa2*sa2); cap = ca1*ca2 - sa1*sa2; /* c+ = cc - ss. */ cam = ca1*ca2 + sa1*sa2; /* c- = cc + ss. */ sap = sa1*ca2 + ca1*sa2; /* s+ = sc + cs. */ sam = sa1*ca2 - ca1*sa2; /* s- = sc - cs. */ r = 0.5*sam*sam*(cam*cam + cap*cap) / (sap*sap*cam*cam); /* rearranged for speed. */ } } return(r); } /******** END SUBROUTINE **********/ /*********************************************************** * the boundary is the face of some voxel * find the the photon's hitting position on the nearest face of the voxel and update the step size. ****/ double FindVoxelFace(double x1, double y1, double z1, double x2, double y2, double z2, double dx, double dy, double dz, double ux, double uy, double uz) { double x_1 = x1 / dx; double y_1 = y1 / dy; double z_1 = z1 / dz; double x_2 = x2 / dx; double y_2 = y2 / dy; double z_2 = z2 / dz; double fx_1 = floor(x_1); double fy_1 = floor(y_1); double fz_1 = floor(z_1); double fx_2 = floor(x_2); double fy_2 = floor(y_2); double fz_2 = floor(z_2); double x = 0, y = 0, z = 0, x0 = 0, y0 = 0, z0 = 0, s = 0; if ((fx_1 != fx_2) && (fy_1 != fy_2) && (fz_1 != fz_2)) { //#10 fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added x = (max2(fx_1, fx_2) - x_1) / ux; y = (max2(fy_1, fy_2) - y_1) / uy; z = (max2(fz_1, fz_2) - z_1) / uz; if (x == min3(x, y, z)) { x0 = max2(fx_1, fx_2); y0 = (x0 - x_1) / ux*uy + y_1; z0 = (x0 - x_1) / ux*uz + z_1; } else if (y == min3(x, y, z)) { y0 = max2(fy_1, fy_2); x0 = (y0 - y_1) / uy*ux + x_1; z0 = (y0 - y_1) / uy*uz + z_1; } else { z0 = max2(fz_1, fz_2); y0 = (z0 - z_1) / uz*uy + y_1; x0 = (z0 - z_1) / uz*ux + x_1; } } else if ((fx_1 != fx_2) && (fy_1 != fy_2)) { //#2 fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added x = (max2(fx_1, fx_2) - x_1) / ux; y = (max2(fy_1, fy_2) - y_1) / uy; if (x == min2(x, y)) { x0 = max2(fx_1, fx_2); y0 = (x0 - x_1) / ux*uy + y_1; z0 = (x0 - x_1) / ux*uz + z_1; } else { y0 = max2(fy_1, fy_2); x0 = (y0 - y_1) / uy*ux + x_1; z0 = (y0 - y_1) / uy*uz + z_1; } } else if ((fy_1 != fy_2) && (fz_1 != fz_2)) { //#3 fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added y = (max2(fy_1, fy_2) - y_1) / uy; z = (max2(fz_1, fz_2) - z_1) / uz; if (y == min2(y, z)) { y0 = max2(fy_1, fy_2); x0 = (y0 - y_1) / uy*ux + x_1; z0 = (y0 - y_1) / uy*uz + z_1; } else { z0 = max2(fz_1, fz_2); x0 = (z0 - z_1) / uz*ux + x_1; y0 = (z0 - z_1) / uz*uy + y_1; } } else if ((fx_1 != fx_2) && (fz_1 != fz_2)) { //#4 fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added x = (max2(fx_1, fx_2) - x_1) / ux; z = (max2(fz_1, fz_2) - z_1) / uz; if (x == min2(x, z)) { x0 = max2(fx_1, fx_2); y0 = (x0 - x_1) / ux*uy + y_1; z0 = (x0 - x_1) / ux*uz + z_1; } else { z0 = max2(fz_1, fz_2); x0 = (z0 - z_1) / uz*ux + x_1; y0 = (z0 - z_1) / uz*uy + y_1; } } else if (fx_1 != fx_2) { //#5 fx_2 = fx_1 + SIGN(fx_2 - fx_1);//added x0 = max2(fx_1, fx_2); y0 = (x0 - x_1) / ux*uy + y_1; z0 = (x0 - x_1) / ux*uz + z_1; } else if (fy_1 != fy_2) { //#6 fy_2 = fy_1 + SIGN(fy_2 - fy_1);//added y0 = max2(fy_1, fy_2); x0 = (y0 - y_1) / uy*ux + x_1; z0 = (y0 - y_1) / uy*uz + z_1; } else { //#7 z0 = max2(fz_1, fz_2); fz_2 = fz_1 + SIGN(fz_2 - fz_1);//added x0 = (z0 - z_1) / uz*ux + x_1; y0 = (z0 - z_1) / uz*uy + y_1; } //s = ( (x0-fx_1)*dx + (y0-fy_1)*dy + (z0-fz_1)*dz )/3; //s = sqrt( SQR((x0-x_1)*dx) + SQR((y0-y_1)*dy) + SQR((z0-z_1)*dz) ); //s = sqrt(SQR(x0-x_1)*SQR(dx) + SQR(y0-y_1)*SQR(dy) + SQR(z0-z_1)*SQR(dz)); s = sqrt(SQR((x0 - x_1)*dx) + SQR((y0 - y_1)*dy) + SQR((z0 - z_1)*dz)); return (s); } }; /* DECLARE FUNCTIONS */ void runMonteCarlo(InputValues* inputValues, ReturnValues* returnValue, int thread, mitk::pa::MonteCarloThreadHandler::Pointer threadHandler); int detector_x = -1; int detector_z = -1; bool interpretAsTime = true; bool simulatePVFC = false; int requestedNumberOfPhotons = 100000; float requestedSimulationTime = 0; // in minutes int concurentThreadsSupported = -1; float yOffset = 0; // in mm bool saveLegacy = false; std::string normalizationFilename; std::string inputFilename; std::string outputFilename; mitk::pa::Probe::Pointer m_PhotoacousticProbe; int main(int argc, char * argv[]) { mitkCommandLineParser parser; // set general information parser.setCategory("MITK-Photoacoustics"); parser.setTitle("Mitk MCxyz"); parser.setDescription("Runs Monte Carlo simulations on inputed tissues."); parser.setContributor("CAI, DKFZ based on code by Jacques and Li"); // 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 I/O parameters"); parser.addArgument( "input", "i", mitkCommandLineParser::InputFile, "Input tissue file", "input tissue file (*.nrrd)", us::Any(), false); parser.addArgument( "output", "o", mitkCommandLineParser::OutputFile, "Output fluence file", "where to save the simulated fluence (*.nrrd)", us::Any(), false); parser.endGroup(); parser.beginGroup("Optional parameters"); parser.addArgument( "verbose", "v", mitkCommandLineParser::Bool, "Verbose Output", "Whether to produce verbose, or rather debug output"); parser.addArgument( "detector-x", "dx", mitkCommandLineParser::Int, "Detector voxel x position", "Determines the x position of the detector voxel (default: -1 = dont use detector voxel)", -1); parser.addArgument( "detector-z", "dz", mitkCommandLineParser::Int, "Detector voxel z position", "Determines the z position of the detector voxel (default: -1 = dont use detector voxel)", -1); parser.addArgument( "number-of-photons", "n", mitkCommandLineParser::Int, "Number of photons", "Specifies the number of photons (default: 100000). Simulation stops after that number. Use -t --timer to define a timer instead"); parser.addArgument( "timer", "t", mitkCommandLineParser::Float, "Simulation time in min", "Specifies the amount of time for simutation (default: 0). Simulation stops after that number of minutes. -n --number-of-photons is the override and default behavior and defines the maximum number of photons instead. If no simulation time or number of photons is specified the file time is taken."); parser.addArgument( "y-offset", "yo", mitkCommandLineParser::Float, "Probe Y-Offset in mm", "Specifies an offset of the photoacoustic probe in the y direction depending on the initial probe position (default: 0) in mm."); parser.addArgument( "jobs", "j", mitkCommandLineParser::Int, "Number of jobs", "Specifies the number of jobs for simutation (default: -1 which starts as many jobs as supported)."); parser.addArgument( "probe-xml", "p", mitkCommandLineParser::InputFile, "Xml definition of the probe", "Specifies the absolute path of the location of the xml definition file of the probe design."); parser.addArgument("normalization-file", "nf", mitkCommandLineParser::InputFile, "Input normalization file", "The input normalization file is used for normalization of the number of photons in the PVFC calculations."); parser.endGroup(); // parse arguments, this method returns a mapping of long argument names and their values std::map parsedArgs = parser.parseArguments(argc, argv); if (parsedArgs.size() == 0) return EXIT_FAILURE; // parse, cast and set required arguments inputFilename = us::any_cast(parsedArgs["input"]); // strip ending inputFilename = inputFilename.substr(0, inputFilename.find("_H.mci")); inputFilename = inputFilename.substr(0, inputFilename.find("_T.bin")); outputFilename = us::any_cast(parsedArgs["output"]); // add .nrrd if not there std::string suffix = ".nrrd"; if (outputFilename.compare(outputFilename.size() - suffix.size(), suffix.size(), suffix) != 0) outputFilename = outputFilename + suffix; // default values for optional arguments // parse, cast and set optional arguments if given if (parsedArgs.count("verbose")) { verbose = us::any_cast(parsedArgs["verbose"]); } if (parsedArgs.count("detector-x")) { detector_x = us::any_cast(parsedArgs["detector-x"]); } if (parsedArgs.count("detector-z")) { detector_z = us::any_cast(parsedArgs["detector-z"]); } if (parsedArgs.count("timer")) { requestedSimulationTime = us::any_cast(parsedArgs["timer"]); if (requestedSimulationTime > 0) interpretAsTime = true; } if (parsedArgs.count("y-offset")) { yOffset = us::any_cast(parsedArgs["y-offset"]); } if (parsedArgs.count("number-of-photons")) { requestedNumberOfPhotons = us::any_cast(parsedArgs["number-of-photons"]); if (requestedNumberOfPhotons > 0) interpretAsTime = false; } if (parsedArgs.count("jobs")) { concurentThreadsSupported = us::any_cast(parsedArgs["jobs"]); } if (parsedArgs.count("probe-xml")) { std::string inputXmlProbeDesign = us::any_cast(parsedArgs["probe-xml"]); m_PhotoacousticProbe = mitk::pa::Probe::New(inputXmlProbeDesign, verbose); if (!m_PhotoacousticProbe->IsValid()) { std::cerr << "Xml File was not valid. Simulation failed." << std::endl; return EXIT_FAILURE; } } if (parsedArgs.count("normalization-file")) { normalizationFilename = us::any_cast(parsedArgs["normalization-file"]); } if (concurentThreadsSupported == 0 || concurentThreadsSupported == -1) { concurentThreadsSupported = std::thread::hardware_concurrency(); if (concurentThreadsSupported == 0) { std::cout << "Could not determine number of available threads. Launching only one." << std::endl; concurentThreadsSupported = 1; } } if (detector_x != -1 && detector_z != -1) { if (verbose) std::cout << "Performing PVFC calculation for x=" << detector_x << " and z=" << detector_z << std::endl; simulatePVFC = true; } else { if (verbose) std::cout << "Will not perform PVFC calculation due to x=" << detector_x << " and/or z=" << detector_z << std::endl; } InputValues allInput = InputValues(); allInput.LoadValues(inputFilename, yOffset, normalizationFilename, simulatePVFC); std::vector allValues(concurentThreadsSupported); std::thread* threads = new std::thread[concurentThreadsSupported]; for (long i = 0; i < concurentThreadsSupported; i++) { ReturnValues* tmp = new ReturnValues(); allValues.push_back(*tmp); } if (verbose) std::cout << "Initializing MonteCarloThreadHandler" << std::endl; long timeMetric; if (interpretAsTime) { if (requestedSimulationTime < mitk::eps) requestedSimulationTime = allInput.simulationTimeFromFile; timeMetric = requestedSimulationTime * 60 * 1000; } else { timeMetric = requestedNumberOfPhotons; } mitk::pa::MonteCarloThreadHandler::Pointer threadHandler = mitk::pa::MonteCarloThreadHandler::New(timeMetric, interpretAsTime); if (simulatePVFC) threadHandler->SetPackageSize(1000); if (verbose) std::cout << "\nStarting simulation ...\n" << std::endl; auto simulationStartTime = std::chrono::system_clock::now(); for (int i = 0; i < concurentThreadsSupported; i++) { threads[i] = std::thread(runMonteCarlo, &allInput, &allValues[i], (i + 1), threadHandler); } for (int i = 0; i < concurentThreadsSupported; i++) { threads[i].join(); } auto simulationFinishTime = std::chrono::system_clock::now(); auto simulationTimeElapsed = simulationFinishTime - simulationStartTime; if (verbose) std::cout << "\n\nFinished simulation\n\n" << std::endl; std::cout << "total time for simulation: " << (int)std::chrono::duration_cast(simulationTimeElapsed).count() << "sec " << std::endl; /**** SAVE Convert data to relative fluence rate [cm^-2] and save. *****/ if (!simulatePVFC) { if (verbose) std::cout << "Allocating memory for normal simulation result ... "; double* finalTotalFluence = (double *)malloc(allInput.totalNumberOfVoxels * sizeof(double)); if (verbose) std::cout << "[OK]" << std::endl; if (verbose) std::cout << "Cleaning memory for normal simulation result ..."; for (int i = 0; i < allInput.totalNumberOfVoxels; i++) { finalTotalFluence[i] = 0; } if (verbose) std::cout << "[OK]" << std::endl; if (verbose) std::cout << "Calculating resulting fluence ... "; double tdx = 0, tdy = 0, tdz = 0; long long tNphotons = 0; for (int t = 0; t < concurentThreadsSupported; t++) { tdx = allInput.xSpacing; tdy = allInput.ySpacing; tdz = allInput.zSpacing; tNphotons += allValues[t].Nphotons; for (int voxelNumber = 0; voxelNumber < allInput.totalNumberOfVoxels; voxelNumber++) { finalTotalFluence[voxelNumber] += allValues[t].totalFluence[voxelNumber]; } } if (verbose) std::cout << "[OK]" << std::endl; std::cout << "total number of photons simulated: " << tNphotons << std::endl; // Normalize deposition (A) to yield fluence rate (F). double temp = tdx*tdy*tdz*tNphotons; for (int i = 0; i < allInput.totalNumberOfVoxels; i++) { finalTotalFluence[i] /= temp*allInput.muaVector[i]; } if (verbose) std::cout << "Saving normal simulated fluence result to " << outputFilename << " ... "; mitk::Image::Pointer resultImage = mitk::Image::New(); mitk::PixelType TPixel = mitk::MakeScalarPixelType(); unsigned int* dimensionsOfImage = new unsigned int[3]; // Copy dimensions dimensionsOfImage[0] = allInput.Ny; dimensionsOfImage[1] = allInput.Nx; dimensionsOfImage[2] = allInput.Nz; resultImage->Initialize(TPixel, 3, dimensionsOfImage); mitk::Vector3D spacing; spacing[0] = allInput.ySpacing; spacing[1] = allInput.xSpacing; spacing[2] = allInput.zSpacing; resultImage->SetSpacing(spacing); resultImage->SetImportVolume(finalTotalFluence, 0, 0, mitk::Image::CopyMemory); resultImage->GetPropertyList()->SetFloatProperty("y-offset", yOffset); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("y-offset")); mitk::IOUtil::Save(resultImage, outputFilename); if (verbose) std::cout << "[OK]" << std::endl; if (verbose) { std::cout << "x spacing = " << tdx << std::endl; std::cout << "y spacing = " << tdy << std::endl; std::cout << "z spacing = " << tdz << std::endl; std::cout << "total number of voxels = " << allInput.totalNumberOfVoxels << std::endl; std::cout << "number of photons = " << (int)tNphotons << std::endl; } } else // if simulate PVFC { if (verbose) std::cout << "Allocating memory for PVFC simulation result ... "; double* detectorFluence = ((double*)malloc(allInput.totalNumberOfVoxels * sizeof(double))); if (verbose) std::cout << "[OK]" << std::endl; if (verbose) std::cout << "Cleaning memory for PVFC simulation result ..."; for (int i = 0; i < allInput.totalNumberOfVoxels; i++) { detectorFluence[i] = 0; } if (verbose) std::cout << "[OK]" << std::endl; if (verbose) std::cout << "Calculating resulting PVFC fluence ... "; double tdx = 0, tdy = 0, tdz = 0; long long tNphotons = 0; long pvfcPhotons = 0; for (int t = 0; t < concurentThreadsSupported; t++) { tdx = allInput.xSpacing; tdy = allInput.ySpacing; tdz = allInput.zSpacing; tNphotons += allValues[t].Nphotons; pvfcPhotons += allValues[t].detectorVoxel->m_NumberPhotonsCurrent; for (int voxelNumber = 0; voxelNumber < allInput.totalNumberOfVoxels; voxelNumber++) { detectorFluence[voxelNumber] += allValues[t].detectorVoxel->fluenceContribution[voxelNumber]; } } if (verbose) std::cout << "[OK]" << std::endl; std::cout << "total number of photons simulated: " << tNphotons << std::endl; // Normalize deposition (A) to yield fluence rate (F). double temp = tdx*tdy*tdz*tNphotons; for (int i = 0; i < allInput.totalNumberOfVoxels; i++) { detectorFluence[i] /= temp*allInput.muaVector[i]; } if (verbose) std::cout << "Saving PVFC ..."; std::stringstream detectorname(""); double detectorX = allValues[0].detectorVoxel->location.x; double detectorY = allValues[0].detectorVoxel->location.y; double detectorZ = allValues[0].detectorVoxel->location.z; detectorname << detectorX << "," << detectorY << "," << detectorZ << "FluenceContribution.nrrd"; // Save the binary file std::string outputFileBase = outputFilename.substr(0, outputFilename.find(".nrrd")); outputFilename = outputFileBase + "_p" + detectorname.str().c_str(); mitk::Image::Pointer pvfcImage = mitk::Image::New(); unsigned int* dimensionsOfPvfcImage = new unsigned int[3]; // Copy dimensions dimensionsOfPvfcImage[0] = allInput.Ny; dimensionsOfPvfcImage[1] = allInput.Nx; dimensionsOfPvfcImage[2] = allInput.Nz; pvfcImage->Initialize(mitk::MakeScalarPixelType(), 3, dimensionsOfPvfcImage); mitk::Vector3D pvfcSpacing; pvfcSpacing[0] = allInput.ySpacing; pvfcSpacing[1] = allInput.xSpacing; pvfcSpacing[2] = allInput.zSpacing; pvfcImage->SetSpacing(pvfcSpacing); pvfcImage->SetImportVolume(detectorFluence, 0, 0, mitk::Image::CopyMemory); pvfcImage->GetPropertyList()->SetFloatProperty("detector-x", detectorX); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("detector-x")); pvfcImage->GetPropertyList()->SetFloatProperty("detector-y", detectorY); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("detector-y")); pvfcImage->GetPropertyList()->SetFloatProperty("detector-z", detectorZ); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("detector-z")); pvfcImage->GetPropertyList()->SetFloatProperty("normalization-factor", allValues[0].detectorVoxel->m_PhotonNormalizationValue); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("normalization-factor")); pvfcImage->GetPropertyList()->SetFloatProperty("simulated-photons", pvfcPhotons); mitk::CoreServices::GetPropertyPersistence()->AddInfo(mitk::PropertyPersistenceInfo::New("simulated-photons")); mitk::IOUtil::Save(pvfcImage, outputFilename); if (verbose) std::cout << "[OK]" << std::endl; if (verbose) { std::cout << "x spacing = " << tdx << std::endl; std::cout << "y spacing = " << tdy << std::endl; std::cout << "z spacing = " << tdz << std::endl; std::cout << "total number of voxels = " << allInput.totalNumberOfVoxels << std::endl; std::cout << "number of photons = " << (int)tNphotons << std::endl; } } exit(EXIT_SUCCESS); } /* end of main */ /* CORE FUNCTION */ void runMonteCarlo(InputValues* inputValues, ReturnValues* returnValue, int thread, mitk::pa::MonteCarloThreadHandler::Pointer threadHandler) { if (verbose) std::cout << "Thread " << thread << ": Locking Mutex ..." << std::endl; if (verbose) std::cout << "[OK]" << std::endl; if (verbose) std::cout << "Initializing ... "; /* Propagation parameters */ double x, y, z; /* photon position */ double ux, uy, uz; /* photon trajectory as cosines */ double uxx, uyy, uzz; /* temporary values used during SPIN */ double s; /* step sizes. s = -log(RND)/mus [cm] */ double sleft; /* dimensionless */ double costheta; /* cos(theta) */ double sintheta; /* sin(theta) */ double cospsi; /* cos(psi) */ double sinpsi; /* sin(psi) */ double psi; /* azimuthal angle */ long photonIterator = 0; /* current photon */ double W; /* photon weight */ double absorb; /* weighted deposited in a step due to absorption */ short photon_status; /* flag = ALIVE=1 or DEAD=0 */ bool sv; /* Are they in the same voxel? */ /* dummy variables */ double rnd; /* assigned random value 0-1 */ double r, phi; /* dummy values */ long i, j; /* dummy indices */ double tempx, tempy, tempz; /* temporary variables, used during photon step. */ int ix, iy, iz; /* Added. Used to track photons */ double temp; /* dummy variable */ int bflag; /* boundary flag: 0 = photon inside volume. 1 = outside volume */ int CNT = 0; returnValue->totalFluence = (double *)malloc(inputValues->totalNumberOfVoxels * sizeof(double)); /* relative fluence rate [W/cm^2/W.delivered] */ if (detector_x != -1 && detector_z != -1) { if (detector_x<0 || detector_x>inputValues->Nx) { std::cout << "Requested detector x position not valid. Needs to be >= 0 and <= " << inputValues->Nx << std::endl; exit(EXIT_FAILURE); } if (detector_z<1 || detector_z>inputValues->Nz) { std::cout << "Requested detector z position not valid. Needs to be > 0 and <= " << inputValues->Nz << std::endl; exit(EXIT_FAILURE); } double photonNormalizationValue = 1 / inputValues->GetNormalizationValue(detector_x, inputValues->Ny / 2, detector_z); returnValue->detectorVoxel = new DetectorVoxel(initLocation(detector_x, inputValues->Ny / 2, detector_z, 0), inputValues->totalNumberOfVoxels, photonNormalizationValue); } /**** ======================== MAJOR CYCLE ============================ *****/ auto duration = std::chrono::system_clock::now().time_since_epoch(); returnValue->RandomGen(0, (std::chrono::duration_cast(duration).count() + thread) % 32000, NULL); /* initiate with seed = 1, or any long integer. */ for (j = 0; j < inputValues->totalNumberOfVoxels; j++) returnValue->totalFluence[j] = 0; // ensure F[] starts empty. /**** RUN Launch N photons, initializing each one before progation. *****/ long photonsToSimulate = 0; do { photonsToSimulate = threadHandler->GetNextWorkPackage(); if (returnValue->detectorVoxel != nullptr) { photonsToSimulate = photonsToSimulate * returnValue->detectorVoxel->m_PhotonNormalizationValue; } if (verbose) MITK_INFO << "Photons to simulate: " << photonsToSimulate; photonIterator = 0L; do { /**** LAUNCH Initialize photon position and trajectory. *****/ photonIterator += 1; /* increment photon count */ W = 1.0; /* set photon weight to one */ photon_status = ALIVE; /* Launch an ALIVE photon */ CNT = 0; /**** SET SOURCE* Launch collimated beam at x,y center.****/ /****************************/ /* Initial position. */ if (m_PhotoacousticProbe.IsNotNull()) { double rnd1 = -1; double rnd2 = -1; double rnd3 = -1; double rnd4 = -1; double rnd5 = -1; double rnd6 = -1; double rnd7 = -1; double rnd8 = -1; while ((rnd1 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); while ((rnd2 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); while ((rnd3 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); while ((rnd4 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); while ((rnd5 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); while ((rnd6 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); while ((rnd7 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); while ((rnd8 = returnValue->RandomGen(1, 0, NULL)) <= 0.0); mitk::pa::LightSource::PhotonInformation info = m_PhotoacousticProbe->GetNextPhoton(rnd1, rnd2, rnd3, rnd4, rnd5, rnd6, rnd7, rnd8); x = info.xPosition; y = yOffset + info.yPosition; z = info.zPosition; ux = info.xAngle; uy = info.yAngle; uz = info.zAngle; if (verbose) std::cout << "Created photon at position (" << x << "|" << y << "|" << z << ") with angles (" << ux << "|" << uy << "|" << uz << ")." << std::endl; } else { /* trajectory */ if (inputValues->launchflag == 1) // manually set launch { x = inputValues->xs; y = inputValues->ys; z = inputValues->zs; ux = inputValues->ux0; uy = inputValues->uy0; uz = inputValues->uz0; } else // use mcflag { if (inputValues->mcflag == 0) // uniform beam { // set launch point and width of beam while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); // avoids rnd = 0 r = inputValues->radius*sqrt(rnd); // radius of beam at launch point while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); // avoids rnd = 0 phi = rnd*2.0*PI; x = inputValues->xs + r*cos(phi); y = inputValues->ys + r*sin(phi); z = inputValues->zs; // set trajectory toward focus while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); // avoids rnd = 0 r = inputValues->waist*sqrt(rnd); // radius of beam at focus while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); // avoids rnd = 0 phi = rnd*2.0*PI; // !!!!!!!!!!!!!!!!!!!!!!! setting input values will braek inputValues->xfocus = r*cos(phi); inputValues->yfocus = r*sin(phi); temp = sqrt((x - inputValues->xfocus)*(x - inputValues->xfocus) + (y - inputValues->yfocus)*(y - inputValues->yfocus) + inputValues->zfocus*inputValues->zfocus); ux = -(x - inputValues->xfocus) / temp; uy = -(y - inputValues->yfocus) / temp; uz = sqrt(1 - ux*ux + uy*uy); } else if (inputValues->mcflag == 5) // Multispectral DKFZ prototype { // set launch point and width of beam while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); //offset in x direction in cm (random) x = (rnd*2.5) - 1.25; while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); double b = ((rnd)-0.5); y = (b > 0 ? yOffset + 1.5 : yOffset - 1.5); z = 0.1; ux = 0; while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); //Angle of beam in y direction uy = sin((rnd*0.42) - 0.21 + (b < 0 ? 1.0 : -1.0) * 0.436); while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); // angle of beam in x direction ux = sin((rnd*0.42) - 0.21); uz = sqrt(1 - ux*ux - uy*uy); } else if (inputValues->mcflag == 4) // Monospectral prototype DKFZ { // set launch point and width of beam while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); //offset in x direction in cm (random) x = (rnd*2.5) - 1.25; while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); double b = ((rnd)-0.5); y = (b > 0 ? yOffset + 0.83 : yOffset - 0.83); z = 0.1; ux = 0; while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); //Angle of beam in y direction uy = sin((rnd*0.42) - 0.21 + (b < 0 ? 1.0 : -1.0) * 0.375); while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); // angle of beam in x direction ux = sin((rnd*0.42) - 0.21); uz = sqrt(1 - ux*ux - uy*uy); } else { // isotropic pt source costheta = 1.0 - 2.0 * returnValue->RandomGen(1, 0, NULL); sintheta = sqrt(1.0 - costheta*costheta); psi = 2.0 * PI * returnValue->RandomGen(1, 0, NULL); cospsi = cos(psi); if (psi < PI) sinpsi = sqrt(1.0 - cospsi*cospsi); else sinpsi = -sqrt(1.0 - cospsi*cospsi); x = inputValues->xs; y = inputValues->ys; z = inputValues->zs; ux = sintheta*cospsi; uy = sintheta*sinpsi; uz = costheta; } } // end use mcflag } /****************************/ /* Get tissue voxel properties of launchpoint. * If photon beyond outer edge of defined voxels, * the tissue equals properties of outermost voxels. * Therefore, set outermost voxels to infinite background value. */ ix = (int)(inputValues->Nx / 2 + x / inputValues->xSpacing); iy = (int)(inputValues->Ny / 2 + y / inputValues->ySpacing); iz = (int)(z / inputValues->zSpacing); if (ix >= inputValues->Nx) ix = inputValues->Nx - 1; if (iy >= inputValues->Ny) iy = inputValues->Ny - 1; if (iz >= inputValues->Nz) iz = inputValues->Nz - 1; if (ix < 0) ix = 0; if (iy < 0) iy = 0; if (iz < 0) iz = 0; /* Get the tissue type of located voxel */ i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy); bflag = 1; // initialize as 1 = inside volume, but later check as photon propagates. if (returnValue->detectorVoxel != nullptr) returnValue->detectorVoxel->recordedPhotonRoute->clear(); /* HOP_DROP_SPIN_CHECK Propagate one photon until it dies as determined by ROULETTE. *******/ do { /**** HOP Take step to new position s = dimensionless stepsize x, uy, uz are cosines of current photon trajectory *****/ while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); /* yields 0 < rnd <= 1 */ sleft = -log(rnd); /* dimensionless step */ CNT += 1; do { // while sleft>0 s = sleft / inputValues->musVector[i]; /* Step size [cm].*/ tempx = x + s*ux; /* Update positions. [cm] */ tempy = y + s*uy; tempz = z + s*uz; sv = returnValue->SameVoxel(x, y, z, tempx, tempy, tempz, inputValues->xSpacing, inputValues->ySpacing, inputValues->zSpacing); if (sv) /* photon in same voxel */ { x = tempx; /* Update positions. */ y = tempy; z = tempz; /**** DROP Drop photon weight (W) into local bin. *****/ absorb = W*(1 - exp(-inputValues->muaVector[i] * s)); /* photon weight absorbed at this step */ W -= absorb; /* decrement WEIGHT by amount absorbed */ // If photon within volume of heterogeneity, deposit energy in F[]. // Normalize F[] later, when save output. if (bflag) { i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy); returnValue->totalFluence[i] += absorb; // only save data if blag==1, i.e., photon inside simulation cube //For each detectorvoxel if (returnValue->detectorVoxel != nullptr) { //Add photon position to the recorded photon route returnValue->detectorVoxel->recordedPhotonRoute->push_back(initLocation(ix, iy, iz, absorb)); //If the photon is currently at the detector position if ((returnValue->detectorVoxel->location.x == ix) && ((returnValue->detectorVoxel->location.y == iy) || (returnValue->detectorVoxel->location.y - 1 == iy)) && (returnValue->detectorVoxel->location.z == iz)) { //For each voxel in the recorded photon route for (unsigned int routeIndex = 0; routeIndex < returnValue->detectorVoxel->recordedPhotonRoute->size(); routeIndex++) { //increment the fluence contribution at that particular position i = (long)(returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).z*inputValues->Ny*inputValues->Nx + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).x*inputValues->Ny + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).y); returnValue->detectorVoxel->fluenceContribution[i] += returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).absorb; } //Clear the recorded photon route returnValue->detectorVoxel->m_NumberPhotonsCurrent++; returnValue->detectorVoxel->recordedPhotonRoute->clear(); } } } /* Update sleft */ sleft = 0; /* dimensionless step remaining */ } else /* photon has crossed voxel boundary */ { /* step to voxel face + "littlest step" so just inside new voxel. */ s = ls + returnValue->FindVoxelFace2(x, y, z, tempx, tempy, tempz, inputValues->xSpacing, inputValues->ySpacing, inputValues->zSpacing, ux, uy, uz); /**** DROP Drop photon weight (W) into local bin. *****/ absorb = W*(1 - exp(-inputValues->muaVector[i] * s)); /* photon weight absorbed at this step */ W -= absorb; /* decrement WEIGHT by amount absorbed */ // If photon within volume of heterogeneity, deposit energy in F[]. // Normalize F[] later, when save output. if (bflag) { // only save data if bflag==1, i.e., photon inside simulation cube //For each detectorvoxel if (returnValue->detectorVoxel != nullptr) { //Add photon position to the recorded photon route returnValue->detectorVoxel->recordedPhotonRoute->push_back(initLocation(ix, iy, iz, absorb)); //If the photon is currently at the detector position if ((returnValue->detectorVoxel->location.x == ix) && ((returnValue->detectorVoxel->location.y == iy) || (returnValue->detectorVoxel->location.y - 1 == iy)) && (returnValue->detectorVoxel->location.z == iz)) { //For each voxel in the recorded photon route for (unsigned int routeIndex = 0; routeIndex < returnValue->detectorVoxel->recordedPhotonRoute->size(); routeIndex++) { //increment the fluence contribution at that particular position i = (long)(returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).z*inputValues->Ny*inputValues->Nx + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).x*inputValues->Ny + returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).y); returnValue->detectorVoxel->fluenceContribution[i] += returnValue->detectorVoxel->recordedPhotonRoute->at(routeIndex).absorb; } //Clear the recorded photon route returnValue->detectorVoxel->m_NumberPhotonsCurrent++; returnValue->detectorVoxel->recordedPhotonRoute->clear(); } } i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy); returnValue->totalFluence[i] += absorb; } /* Update sleft */ sleft -= s*inputValues->musVector[i]; /* dimensionless step remaining */ if (sleft <= ls) sleft = 0; /* Update positions. */ x += s*ux; y += s*uy; z += s*uz; // pointers to voxel containing optical properties ix = (int)(inputValues->Nx / 2 + x / inputValues->xSpacing); iy = (int)(inputValues->Ny / 2 + y / inputValues->ySpacing); iz = (int)(z / inputValues->zSpacing); bflag = 1; // Boundary flag. Initialize as 1 = inside volume, then check. if (inputValues->boundaryflag == 0) { // Infinite medium. // Check if photon has wandered outside volume. // If so, set tissue type to boundary value, but let photon wander. // Set blag to zero, so DROP does not deposit energy. if (iz >= inputValues->Nz) { iz = inputValues->Nz - 1; bflag = 0; } if (ix >= inputValues->Nx) { ix = inputValues->Nx - 1; bflag = 0; } if (iy >= inputValues->Ny) { iy = inputValues->Ny - 1; bflag = 0; } if (iz < 0) { iz = 0; bflag = 0; } if (ix < 0) { ix = 0; bflag = 0; } if (iy < 0) { iy = 0; bflag = 0; } } else if (inputValues->boundaryflag == 1) { // Escape at boundaries if (iz >= inputValues->Nz) { iz = inputValues->Nz - 1; photon_status = DEAD; sleft = 0; } if (ix >= inputValues->Nx) { ix = inputValues->Nx - 1; photon_status = DEAD; sleft = 0; } if (iy >= inputValues->Ny) { iy = inputValues->Ny - 1; photon_status = DEAD; sleft = 0; } if (iz < 0) { iz = 0; photon_status = DEAD; sleft = 0; } if (ix < 0) { ix = 0; photon_status = DEAD; sleft = 0; } if (iy < 0) { iy = 0; photon_status = DEAD; sleft = 0; } } else if (inputValues->boundaryflag == 2) { // Escape at top surface, no x,y bottom z boundaries if (iz >= inputValues->Nz) { iz = inputValues->Nz - 1; bflag = 0; } if (ix >= inputValues->Nx) { ix = inputValues->Nx - 1; bflag = 0; } if (iy >= inputValues->Ny) { iy = inputValues->Ny - 1; bflag = 0; } if (iz < 0) { iz = 0; photon_status = DEAD; sleft = 0; } if (ix < 0) { ix = 0; bflag = 0; } if (iy < 0) { iy = 0; bflag = 0; } } // update pointer to tissue type i = (long)(iz*inputValues->Ny*inputValues->Nx + ix*inputValues->Ny + iy); } //(sv) /* same voxel */ } while (sleft > 0); //do...while /**** SPIN Scatter photon into new trajectory defined by theta and psi. Theta is specified by cos(theta), which is determined based on the Henyey-Greenstein scattering function. Convert theta and psi into cosines ux, uy, uz. *****/ /* Sample for costheta */ while ((rnd = returnValue->RandomGen(1, 0, NULL)) <= 0.0); if (inputValues->gVector[i] == 0.0) { costheta = 2.0 * rnd - 1.0; } else { double temp = (1.0 - inputValues->gVector[i] * inputValues->gVector[i]) / (1.0 - inputValues->gVector[i] + 2 * inputValues->gVector[i] * rnd); costheta = (1.0 + inputValues->gVector[i] * inputValues->gVector[i] - temp*temp) / (2.0*inputValues->gVector[i]); } sintheta = sqrt(1.0 - costheta*costheta); /* sqrt() is faster than sin(). */ /* Sample psi. */ psi = 2.0*PI*returnValue->RandomGen(1, 0, NULL); cospsi = cos(psi); if (psi < PI) sinpsi = sqrt(1.0 - cospsi*cospsi); /* sqrt() is faster than sin(). */ else sinpsi = -sqrt(1.0 - cospsi*cospsi); /* New trajectory. */ if (1 - fabs(uz) <= ONE_MINUS_COSZERO) { /* close to perpendicular. */ uxx = sintheta * cospsi; uyy = sintheta * sinpsi; uzz = costheta * SIGN(uz); /* SIGN() is faster than division. */ } else { /* usually use this option */ temp = sqrt(1.0 - uz * uz); uxx = sintheta * (ux * uz * cospsi - uy * sinpsi) / temp + ux * costheta; uyy = sintheta * (uy * uz * cospsi + ux * sinpsi) / temp + uy * costheta; uzz = -sintheta * cospsi * temp + uz * costheta; } /* Update trajectory */ ux = uxx; uy = uyy; uz = uzz; /**** CHECK ROULETTE If photon weight below THRESHOLD, then terminate photon using Roulette technique. Photon has CHANCE probability of having its weight increased by factor of 1/CHANCE, and 1-CHANCE probability of terminating. *****/ if (W < THRESHOLD) { if (returnValue->RandomGen(1, 0, NULL) <= CHANCE) W /= CHANCE; else photon_status = DEAD; } } while (photon_status == ALIVE); /* end STEP_CHECK_HOP_SPIN */ /* if ALIVE, continue propagating */ /* If photon DEAD, then launch new photon. */ } while (photonIterator < photonsToSimulate); /* end RUN */ returnValue->Nphotons += photonsToSimulate; } while (photonsToSimulate > 0); if (verbose) std::cout << "------------------------------------------------------" << std::endl; if (verbose) std::cout << "Thread " << thread << " is finished." << std::endl; } diff --git a/Modules/PhotoacousticsLib/MCxyz/files.cmake b/Modules/PhotoacousticsLib/MitkMCxyz/files.cmake similarity index 50% rename from Modules/PhotoacousticsLib/MCxyz/files.cmake rename to Modules/PhotoacousticsLib/MitkMCxyz/files.cmake index bd0800f8d7..284bae1500 100644 --- a/Modules/PhotoacousticsLib/MCxyz/files.cmake +++ b/Modules/PhotoacousticsLib/MitkMCxyz/files.cmake @@ -1,3 +1,3 @@ set(CPP_FILES - MCxyz.cpp + MitkMCxyz.cpp ) diff --git a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp index cb44c48eb0..034e1aa585 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp @@ -1,317 +1,325 @@ /*=================================================================== 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) { { 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]; for (unsigned int index = 0; index < size; index++) { absorptionArray[index] = parameters->GetBackgroundAbsorption(); scatteringArray[index] = parameters->GetBackgroundScattering(); anisotropyArray[index] = parameters->GetBackgroundAnisotropy(); segmentationArray[index] = SegmentationType::BACKGROUND; } m_AbsorptionVolume = Volume::New(absorptionArray, xDim, yDim, zDim); m_ScatteringVolume = Volume::New(scatteringArray, xDim, yDim, zDim); m_AnisotropyVolume = Volume::New(anisotropyArray, xDim, yDim, zDim); m_SegmentationVolume = Volume::New(segmentationArray, xDim, yDim, zDim); } 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()); if (m_TissueParameters->GetDoVolumeSmoothing()) AddDoubleProperty("sigma", m_TissueParameters->GetVolumeSmoothingSigma()); else AddDoubleProperty("sigma", 0); AddDoubleProperty("standardTissueAbsorption", m_TissueParameters->GetBackgroundAbsorption()); 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; } 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; MITK_INFO << "Initializing image..."; resultImage->Initialize(TPixel, 4, dimensionsOfImage, 1); MITK_INFO << "Initializing image...[Done]"; mitk::Vector3D spacing; spacing.Fill(m_TissueParameters->GetVoxelSpacingInCentimeters()); resultImage->SetSpacing(spacing); MITK_INFO << "Set Import Volumes..."; //Copy memory, deal with cleaning up memory yourself! 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); MITK_INFO << "Set Import Volumes...[Done]"; 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, 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()); } 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; } diff --git a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp index c20729f4e5..01ad8bc6a0 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAVolume.cpp @@ -1,115 +1,115 @@ /*=================================================================== 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 "mitkPAVolume.h" #include #include #include mitk::pa::Volume::Volume(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim) { if (data == nullptr) mitkThrow() << "You may not initialize a mitk::Volume with a nullptr"; m_InternalMitkImage = mitk::Image::New(); unsigned int* dimensions = new unsigned int[NUMBER_OF_SPATIAL_DIMENSIONS]; dimensions[0] = yDim; dimensions[1] = xDim; dimensions[2] = zDim; m_XDim = xDim; m_YDim = yDim; m_ZDim = zDim; mitk::PixelType pixelType = mitk::MakeScalarPixelType(); m_InternalMitkImage->Initialize(pixelType, NUMBER_OF_SPATIAL_DIMENSIONS, dimensions); - m_InternalMitkImage->SetImportVolume(data, Image::ImportMemoryManagementType::RtlCopyMemory); + m_InternalMitkImage->SetImportVolume(data, Image::ImportMemoryManagementType::CopyMemory); m_FastAccessDataPointer = GetData(); delete data; } mitk::pa::Volume::~Volume() { m_InternalMitkImage = nullptr; } mitk::pa::Volume::Pointer mitk::pa::Volume::New(double* data, unsigned int xDim, unsigned int yDim, unsigned int zDim) { mitk::pa::Volume::Pointer smartPtr = new mitk::pa::Volume(data, xDim, yDim, zDim); smartPtr->UnRegister(); return smartPtr; } mitk::Image::Pointer mitk::pa::Volume::AsMitkImage() { return m_InternalMitkImage; } mitk::pa::Volume::Pointer mitk::pa::Volume::DeepCopy() { long length = GetXDim()*GetYDim()*GetZDim(); double* data = new double[length]; memcpy(data, GetData(), length * sizeof(double)); return mitk::pa::Volume::New(data, GetXDim(), GetYDim(), GetZDim()); } double mitk::pa::Volume::GetData(unsigned int x, unsigned int y, unsigned int z) { return m_FastAccessDataPointer[GetIndex(x, y, z)]; } void mitk::pa::Volume::SetData(double data, unsigned int x, unsigned int y, unsigned int z) { m_FastAccessDataPointer[GetIndex(x, y, z)] = data; } unsigned int mitk::pa::Volume::GetXDim() { return m_XDim; } unsigned int mitk::pa::Volume::GetYDim() { return m_YDim; } unsigned int mitk::pa::Volume::GetZDim() { return m_ZDim; } double* mitk::pa::Volume::GetData() const { mitk::ImageWriteAccessor imgRead(m_InternalMitkImage, m_InternalMitkImage->GetVolumeData()); return (double*)imgRead.GetData(); } int mitk::pa::Volume::GetIndex(unsigned int x, unsigned int y, unsigned int z) { #ifdef _DEBUG if (x < 0 || x >(GetXDim() - 1) || y < 0 || y >(GetYDim() - 1) || z < 0 || z >(GetZDim() - 1)) { MITK_ERROR << "Index out of bounds at " << x << "|" << y << "|" << z; mitkThrow() << "Index out of bounds exception!"; } #endif return z * m_XDim * m_YDim + x * m_YDim + y; } diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPANoiseGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPANoiseGenerator.cpp index d29eaa90c7..1bbf001318 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPANoiseGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPANoiseGenerator.cpp @@ -1,68 +1,67 @@ /*=================================================================== 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 "mitkPANoiseGenerator.h" #include "mitkPAIOUtil.h" #include "mitkPAExceptions.h" #include #include #include void mitk::pa::NoiseGenerator::ApplyNoiseModel(mitk::pa::Volume::Pointer image, double detectorNoise, double speckleNoise) { if (detectorNoise < 0 || speckleNoise < 0) throw mitk::pa::InvalidInputException("detectorNoise must be >= 0 and speckleNoise must be >= 0"); if (detectorNoise < mitk::eps && speckleNoise < mitk::eps) return; std::mt19937 rng; std::random_device randomDevice; if (randomDevice.entropy() > mitk::eps) { rng.seed(randomDevice()); } else { rng.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count()); } std::normal_distribution<> detector(detectorNoise / 2, detectorNoise / 2); std::normal_distribution<> speckle(1, speckleNoise); - int negativecounter = 0; - int length = image->GetXDim() * image->GetYDim() * image->GetZDim(); + unsigned int negativecounter = 0; double* data = image->GetData(); - for (int x = 0, xLength = image->GetXDim(); x < xLength; x++) - for (int y = 0, yLength = image->GetYDim(); y < yLength; y++) - for (int z = 0, zLength = image->GetZDim(); z < zLength; z++) + for (unsigned int x = 0, xLength = image->GetXDim(); x < xLength; x++) + for (unsigned int y = 0, yLength = image->GetYDim(); y < yLength; y++) + for (unsigned int z = 0, zLength = image->GetZDim(); z < zLength; z++) { double additiveNoise = detector(rng); double multiplicativeNoise = speckle(rng); double newSignal = (data[image->GetIndex(x, y, z)] + additiveNoise)*multiplicativeNoise; if (newSignal <= mitk::eps) { newSignal = mitk::eps; negativecounter++; } data[image->GetIndex(x, y, z)] = newSignal; } } diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPASimulationBatchGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPASimulationBatchGenerator.cpp index c4c5886ee0..521d2d2fbf 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPASimulationBatchGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPASimulationBatchGenerator.cpp @@ -1,92 +1,91 @@ /*=================================================================== 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 "mitkPASimulationBatchGenerator.h" #include #include #ifdef __linux__ #include #include #else #include #endif mitk::pa::SimulationBatchGenerator::SimulationBatchGenerator() { } mitk::pa::SimulationBatchGenerator::~SimulationBatchGenerator() { } std::string mitk::pa::SimulationBatchGenerator::GetVolumeNumber( SimulationBatchGeneratorParameters::Pointer parameters) { std::string volumeNumber = std::to_string(parameters->GetVolumeIndex()); volumeNumber = std::string((3 - volumeNumber.length()), '0') + volumeNumber; return volumeNumber; } std::string mitk::pa::SimulationBatchGenerator::GetOutputFolderName( SimulationBatchGeneratorParameters::Pointer parameters) { return (parameters->GetTissueName() + GetVolumeNumber(parameters)); } std::string mitk::pa::SimulationBatchGenerator::CreateBatchSimulationString( SimulationBatchGeneratorParameters::Pointer parameters) { std::string outputFolderName = GetOutputFolderName(parameters); std::string savePath = outputFolderName + ".nrrd"; std::stringstream batchstring; - double d = 0; for (double d = parameters->GetYOffsetLowerThresholdInCentimeters(); d <= parameters->GetYOffsetUpperThresholdInCentimeters() + 1e-5; d += parameters->GetYOffsetStepInCentimeters()) { batchstring << parameters->GetBinaryPath() << " -p PROBE_DESIGN.xml -i " << savePath << " -o " << outputFolderName << "/" << parameters->GetTissueName() << GetVolumeNumber(parameters) << "_yo" << round(d * 100) / 100 << ".nrrd" << " -yo " << round(d * 100) / 100 << " -n " << parameters->GetNumberOfPhotons() << "\n"; } return batchstring.str(); } void mitk::pa::SimulationBatchGenerator::WriteBatchFileAndSaveTissueVolume( SimulationBatchGeneratorParameters::Pointer parameters, mitk::Image::Pointer tissueVolume) { std::string outputFolderName = parameters->GetNrrdFilePath() + GetOutputFolderName(parameters); std::string savePath = outputFolderName + ".nrrd"; mitk::IOUtil::Save(tissueVolume, savePath); std::string filenameAllSimulation = "simulate_all"; #ifdef __linux__ mkdir(outputFolderName.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); filenameAllSimulation += ".sh"; #else mkdir(outputFolderName.c_str()); filenameAllSimulation += ".bat"; #endif std::ofstream fileAllSimulation(parameters->GetNrrdFilePath() + "/" + filenameAllSimulation, std::ios_base::app); if (fileAllSimulation.is_open()) { fileAllSimulation << CreateBatchSimulationString(parameters); fileAllSimulation.close(); } } diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPASlicedVolumeGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPASlicedVolumeGenerator.cpp index d7e8de8ebf..ce6ea25487 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPASlicedVolumeGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPASlicedVolumeGenerator.cpp @@ -1,112 +1,118 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPASlicedVolumeGenerator.h" #include mitk::pa::SlicedVolumeGenerator::SlicedVolumeGenerator(int centralYSlice, bool precorrect, Volume::Pointer precorrectionVolume, bool inverse) { m_CentralYSlice = centralYSlice; m_Precorrect = precorrect; m_Inverse = inverse; m_PrecorrectionVolume = nullptr; if (m_Precorrect) + { m_PrecorrectionVolume = precorrectionVolume; + } } mitk::pa::SlicedVolumeGenerator::~SlicedVolumeGenerator() { m_CentralYSlice = -1; m_Precorrect = false; m_PrecorrectionVolume = nullptr; } mitk::pa::Volume::Pointer mitk::pa::SlicedVolumeGenerator::GetSlicedFluenceImageFromComposedVolume( ComposedVolume::Pointer composedVolume) { int fluenceComponents = composedVolume->GetNumberOfFluenceComponents(); int xDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetXDim(); int zDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetZDim(); double* imageArray = new double[xDim*zDim*fluenceComponents]; for (int fluenceComponentIdx = 0; fluenceComponentIdx < fluenceComponents; fluenceComponentIdx++) for (int z = 0; z < zDim; z++) for (int x = 0; x < xDim; x++) { int index = z * xDim * fluenceComponents + x * fluenceComponents + fluenceComponentIdx; imageArray[index] = composedVolume->GetFluenceValue(fluenceComponentIdx, x, m_CentralYSlice, z); if (m_Precorrect) + { imageArray[index] = imageArray[index] / m_PrecorrectionVolume->GetData(x, m_CentralYSlice, z); + } if (m_Inverse) + { if (abs(imageArray[index] - 0) >= mitk::eps) imageArray[index] = 1 / imageArray[index]; else imageArray[index] = INFINITY; + } } return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim); } mitk::pa::Volume::Pointer mitk::pa::SlicedVolumeGenerator::GetSlicedSignalImageFromComposedVolume( ComposedVolume::Pointer composedVolume) { int fluenceComponents = composedVolume->GetNumberOfFluenceComponents(); int xDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetXDim(); int zDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetZDim(); double* imageArray = new double[xDim*zDim*fluenceComponents]; for (int fluenceComponentIdx = 0; fluenceComponentIdx < fluenceComponents; fluenceComponentIdx++) for (int z = 0; z < zDim; z++) for (int x = 0; x < xDim; x++) { int y = m_CentralYSlice + composedVolume->GetYOffsetForFluenceComponentInPixels(fluenceComponentIdx); imageArray[z * xDim * fluenceComponents + x * fluenceComponents + fluenceComponentIdx] = composedVolume->GetFluenceValue(fluenceComponentIdx, x, m_CentralYSlice, z) * composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetData(x, y, z); } return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim); } mitk::pa::Volume::Pointer mitk::pa::SlicedVolumeGenerator::GetSlicedGroundTruthImageFromComposedVolume( ComposedVolume::Pointer composedVolume) { int fluenceComponents = composedVolume->GetNumberOfFluenceComponents(); int xDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetXDim(); int zDim = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetZDim(); double* imageArray = new double[xDim*zDim*fluenceComponents]; for (int fluenceComponentIdx = 0; fluenceComponentIdx < fluenceComponents; fluenceComponentIdx++) for (int z = 0; z < zDim; z++) for (int x = 0; x < xDim; x++) { int y = m_CentralYSlice + composedVolume->GetYOffsetForFluenceComponentInPixels(fluenceComponentIdx); imageArray[z * xDim * fluenceComponents + x * fluenceComponents + fluenceComponentIdx] = composedVolume->GetGroundTruthVolume()->GetAbsorptionVolume()->GetData(x, y, z); } return mitk::pa::Volume::New(imageArray, xDim, fluenceComponents, zDim); } diff --git a/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp b/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp index 306f7b9872..6f0b374dfa 100644 --- a/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp +++ b/Modules/PhotoacousticsLib/src/IO/mitkPAIOUtil.cpp @@ -1,250 +1,253 @@ #include "mitkPAIOUtil.h" #include "mitkIOUtil.h" #include "mitkImageReadAccessor.h" -#include "boost/filesystem.hpp" - #include #include #include #include "mitkPAComposedVolume.h" #include "mitkPASlicedVolumeGenerator.h" #include "mitkPANoiseGenerator.h" #include "mitkPAVolumeManipulator.h" #include +#include +#include static std::vector splitString(const std::string &s, const char* delim) { std::vector elems; std::stringstream ss(s); std::string item; while (std::getline(ss, item, *delim)) { int numb; std::stringstream(item) >> numb; elems.push_back(numb); } return elems; } bool mitk::pa::IOUtil::DoesFileHaveEnding(std::string const &fullString, std::string const &ending) { if (fullString.length() == 0 || ending.length() == 0 || fullString.length() < ending.length()) return false; return (0 == fullString.compare(fullString.length() - ending.length(), ending.length(), ending)); } mitk::pa::IOUtil::IOUtil() {} mitk::pa::IOUtil::~IOUtil() {} mitk::pa::Volume::Pointer mitk::pa::IOUtil::LoadNrrd(std::string filename, double blur) { if (filename.empty() || filename == "") return nullptr; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(filename); if (inputImage.IsNull()) return nullptr; int xDim = inputImage->GetDimensions()[1]; int yDim = inputImage->GetDimensions()[0]; int zDim = inputImage->GetDimensions()[2]; mitk::ImageReadAccessor readAccess(inputImage, inputImage->GetVolumeData(0)); double* dataArray = new double[xDim*yDim*zDim]; const double* srcData = (const double*)readAccess.GetData(); memcpy(dataArray, srcData, xDim*yDim*zDim * sizeof(double)); auto returnImage = mitk::pa::Volume::New(dataArray, xDim, yDim, zDim); mitk::pa::VolumeManipulator::GaussianBlur3D(returnImage, blur); return returnImage; } std::map mitk::pa::IOUtil::LoadFluenceContributionMaps(std::string foldername, double blur, int* progress, bool doLog10) { std::map resultMap; - boost::filesystem::directory_iterator end_itr; - for (boost::filesystem::directory_iterator itr(foldername); itr != end_itr; itr++) + itk::Directory::Pointer directoryHandler = itk::Directory::New(); + directoryHandler->Load(foldername.c_str()); + for (unsigned int fileIndex = 0, numFiles = directoryHandler->GetNumberOfFiles(); fileIndex < numFiles; ++fileIndex) { - if (boost::filesystem::is_directory(itr->status())) + std::string filename = std::string(directoryHandler->GetFile(fileIndex)); + if (itksys::SystemTools::FileIsDirectory(filename)) continue; - std::string filename = itr->path().string(); if (!DoesFileHaveEnding(filename, ".nrrd")) continue; size_t s = filename.find("_p"); size_t e = filename.find("Fluence", s); std::string sub = filename.substr(s + 2, e - s - 2); std::vector coords = splitString(sub, ","); if (coords.size() != 3) { MITK_ERROR << "Some of the data to read was corrupted or did not match the " << "naming pattern *_pN,N,NFluence*.nrrd"; mitkThrow() << "Some of the data to read was corrupted or did not match the" << " naming pattern *_pN,N,NFluence*.nrrd"; } else { MITK_DEBUG << "Extracted coords: " << coords[0] << "|" << coords[1] << "|" << coords[2] << " from string " << sub; - Volume::Pointer nrrdFile = LoadNrrd(filename, blur); + Volume::Pointer nrrdFile = LoadNrrd(foldername + filename, blur); if (doLog10) VolumeManipulator::Log10Image(nrrdFile); resultMap[Position{ coords[0], coords[2] }] = nrrdFile; *progress = *progress + 1; } } return resultMap; } int mitk::pa::IOUtil::GetNumberOfNrrdFilesInDirectory(std::string directory) { return GetListOfAllNrrdFilesInDirectory(directory).size(); } std::vector mitk::pa::IOUtil::GetListOfAllNrrdFilesInDirectory(std::string directory, bool keepFileFormat) { std::vector filenames; - boost::filesystem::directory_iterator end_itr; - for (boost::filesystem::directory_iterator itr(directory); itr != end_itr; itr++) + itk::Directory::Pointer directoryHandler = itk::Directory::New(); + directoryHandler->Load(directory.c_str()); + for (unsigned int fileIndex = 0, numFiles = directoryHandler->GetNumberOfFiles(); fileIndex < numFiles; ++fileIndex) { - if (boost::filesystem::is_directory(itr->status())) + std::string filename = std::string(directoryHandler->GetFile(fileIndex)); + if (itksys::SystemTools::FileIsDirectory(filename)) continue; - std::string filename = itr->path().filename().string(); if (!DoesFileHaveEnding(filename, ".nrrd")) continue; if (keepFileFormat) { filenames.push_back(filename); } else { filenames.push_back(filename.substr(0, filename.size() - 5)); } } return filenames; } std::vector mitk::pa::IOUtil::GetAllChildfoldersFromFolder(std::string folderPath) { std::vector returnVector; - boost::filesystem::directory_iterator end_itr; - for (boost::filesystem::directory_iterator itr(folderPath + "/"); itr != end_itr; itr++) + itksys::Directory directoryHandler; + directoryHandler.Load(folderPath.c_str()); + for (unsigned int fileIndex = 2, numFiles = directoryHandler.GetNumberOfFiles(); fileIndex < numFiles; ++fileIndex) { - std::string filename = itr->path().filename().string(); - if (boost::filesystem::is_directory(itr->status())) + std::string filename = folderPath + "/" + std::string(directoryHandler.GetFile(fileIndex)); + if (itksys::SystemTools::FileIsDirectory(filename)) { - returnVector.push_back(folderPath + "/" + filename); + returnVector.push_back(filename); continue; } //If there is a nrrd file in the directory we assume that a bottom level directory was chosen. if (DoesFileHaveEnding(filename, ".nrrd")) { returnVector.clear(); returnVector.push_back(folderPath); return returnVector; } } return returnVector; } mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::IOUtil::LoadInSilicoTissueVolumeFromNrrdFile(std::string nrrdFile) { MITK_INFO << "Initializing ComposedVolume by nrrd..."; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(nrrdFile); auto tissueParameters = TissueGeneratorParameters::New(); unsigned int xDim = inputImage->GetDimensions()[1]; unsigned int yDim = inputImage->GetDimensions()[0]; unsigned int zDim = inputImage->GetDimensions()[2]; tissueParameters->SetXDim(xDim); tissueParameters->SetYDim(yDim); tissueParameters->SetZDim(zDim); double xSpacing = inputImage->GetGeometry(0)->GetSpacing()[1]; double ySpacing = inputImage->GetGeometry(0)->GetSpacing()[0]; double zSpacing = inputImage->GetGeometry(0)->GetSpacing()[2]; if ((xSpacing - ySpacing) > mitk::eps || (xSpacing - zSpacing) > mitk::eps || (ySpacing - zSpacing) > mitk::eps) { throw mitk::Exception("Cannot handle unequal spacing."); } tissueParameters->SetVoxelSpacingInCentimeters(xSpacing); mitk::PropertyList::Pointer propertyList = inputImage->GetPropertyList(); mitk::ImageReadAccessor readAccess0(inputImage, inputImage->GetVolumeData(0)); double* m_AbsorptionArray = new double[xDim*yDim*zDim]; memcpy(m_AbsorptionArray, readAccess0.GetData(), xDim*yDim*zDim * sizeof(double)); auto absorptionVolume = Volume::New(m_AbsorptionArray, xDim, yDim, zDim); mitk::ImageReadAccessor readAccess1(inputImage, inputImage->GetVolumeData(1)); double* m_ScatteringArray = new double[xDim*yDim*zDim]; memcpy(m_ScatteringArray, readAccess1.GetData(), xDim*yDim*zDim * sizeof(double)); auto scatteringVolume = Volume::New(m_ScatteringArray, xDim, yDim, zDim); mitk::ImageReadAccessor readAccess2(inputImage, inputImage->GetVolumeData(2)); double* m_AnisotropyArray = new double[xDim*yDim*zDim]; memcpy(m_AnisotropyArray, readAccess2.GetData(), xDim*yDim*zDim * sizeof(double)); auto anisotropyVolume = Volume::New(m_AnisotropyArray, xDim, yDim, zDim); Volume::Pointer segmentationVolume; if (inputImage->GetDimension() == 4) { mitk::ImageReadAccessor readAccess3(inputImage, inputImage->GetVolumeData(3)); double* m_SegmentationArray = new double[xDim*yDim*zDim]; memcpy(m_SegmentationArray, readAccess3.GetData(), xDim*yDim*zDim * sizeof(double)); segmentationVolume = Volume::New(m_SegmentationArray, xDim, yDim, zDim); } return mitk::pa::InSilicoTissueVolume::New(absorptionVolume, scatteringVolume, anisotropyVolume, segmentationVolume, tissueParameters, propertyList); } mitk::pa::FluenceYOffsetPair::Pointer mitk::pa::IOUtil::LoadFluenceSimulation(std::string fluenceSimulation) { MITK_INFO << "Adding slice..."; mitk::Image::Pointer inputImage = mitk::IOUtil::LoadImage(fluenceSimulation); mitk::ImageReadAccessor readAccess0(inputImage, inputImage->GetVolumeData(0)); unsigned int xDim = inputImage->GetDimensions()[1]; unsigned int yDim = inputImage->GetDimensions()[0]; unsigned int zDim = inputImage->GetDimensions()[2]; int size = xDim*yDim*zDim; double* fluenceArray = new double[size]; memcpy(fluenceArray, readAccess0.GetData(), size * sizeof(double)); auto yOffsetProperty = inputImage->GetProperty("y-offset"); if (yOffsetProperty.IsNull()) mitkThrow() << "No \"y-offset\" property found in fluence file!"; std::string yOff = yOffsetProperty->GetValueAsString(); MITK_INFO << "Reading y Offset: " << yOff; #ifdef __linux__ std::replace(yOff.begin(), yOff.end(), '.', ','); #endif // __linux__ double yOffset = std::stod(yOff); MITK_INFO << "Converted offset " << yOffset; return mitk::pa::FluenceYOffsetPair::New(mitk::pa::Volume::New(fluenceArray, xDim, yDim, zDim), yOffset); } diff --git a/Modules/PhotoacousticsLib/src/Utils/ProbeDesign/mitkPALightSource.cpp b/Modules/PhotoacousticsLib/src/Utils/ProbeDesign/mitkPALightSource.cpp index 437da535a5..6ba871b985 100644 --- a/Modules/PhotoacousticsLib/src/Utils/ProbeDesign/mitkPALightSource.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/ProbeDesign/mitkPALightSource.cpp @@ -1,414 +1,414 @@ /*=================================================================== 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 "mitkPALightSource.h" #include "math.h" mitk::pa::LightSource::LightSource() : m_IsValid(false) { } mitk::pa::LightSource::LightSource(TiXmlElement* element, bool verbose) : m_IsValid(true), m_Verbose(verbose) { ParseEnergy(element); ParsePhotonSpawnArea(element); ParsePhotonDirection(element); if (m_IsValid) { if (m_Verbose) std::cout << "Successfully created LightSource" << std::endl; } else { if (m_Verbose) std::cout << "Failed creating LightSource." << std::endl; } } mitk::pa::LightSource::~LightSource() { } mitk::pa::LightSource::TransformResult mitk::pa::LightSource::BoxMuellerTransform(double u1, double u2, double mu, double sigma) { TransformResult result; result.z0 = sqrt(-2.0 * log(u1)) * cos(TWO_PI * u2) * sigma + mu; result.z1 = sqrt(-2.0 * log(u1)) * sin(TWO_PI * u2) * sigma + mu; return result; } void mitk::pa::LightSource::ParsePhotonDirection(TiXmlElement* element) { TiXmlElement* direction = element->FirstChildElement(XML_TAG_PHOTON_DIRECTION); if (direction) { ParseAngle(direction, XML_TAG_X_ANGLE); } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_X_ANGLE << "\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl; m_AngleXMinimum = 0; m_AngleXMaximum = 0; m_AngleXMode = DistributionMode::UNIFORM; } direction = element->FirstChildElement(XML_TAG_PHOTON_DIRECTION); if (direction) { ParseAngle(direction, XML_TAG_Y_ANGLE); } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_Y_ANGLE << "\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl; m_AngleYMinimum = 0; m_AngleYMaximum = 0; m_AngleYMode = DistributionMode::UNIFORM; } } void mitk::pa::LightSource::ParseAngle(TiXmlElement* direction, std::string angle) { double minimum; double maximum; - DistributionMode mode; + DistributionMode mode = DistributionMode::GAUSSIAN; if (m_Verbose) std::cout << "Parsing " << angle << std::endl; TiXmlElement* angleElement = direction->FirstChildElement(angle); if (angleElement) { TiXmlElement* angleMin = angleElement->FirstChildElement(XML_TAG_MINIMUM); if (angleMin) { std::string angleMinText = angleMin->GetText(); minimum = std::stod(angleMinText); if (m_Verbose) std::cout << "Setting min=" << minimum << std::endl; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_MINIMUM << "\" tag in xml. Setting min=0" << std::endl; minimum = 0; } TiXmlElement* angleMax = angleElement->FirstChildElement(XML_TAG_MAXIMUM); if (angleMax) { std::string angleMaxText = angleMax->GetText(); maximum = std::stod(angleMaxText); if (m_Verbose) std::cout << "Setting max=" << maximum << std::endl; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_MAXIMUM << "\" tag in xml. Setting max=0" << std::endl; maximum = 0; } TiXmlElement* angleMode = angleElement->FirstChildElement(XML_TAG_MODE); if (angleMode) { std::string angleModeText = angleMode->GetText(); if (strcmp("UNIFORM", angleModeText.c_str()) == 0) { mode = DistributionMode::UNIFORM; if (m_Verbose) std::cout << "Setting mode=UNIFORM" << std::endl; } else if (strcmp("GAUSSIAN", angleModeText.c_str()) == 0) { mode = DistributionMode::GAUSSIAN; if (m_Verbose) std::cout << "Setting mode=GAUSSIAN" << std::endl; } } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_MODE << "\" tag in xml. Setting mode=UNIFORM" << std::endl; mode = DistributionMode::UNIFORM; } } else { if (m_Verbose) std::cerr << "No \"" << angle << "\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl; maximum = 0; minimum = 0; mode = DistributionMode::UNIFORM; } if (strcmp(XML_TAG_X_ANGLE.c_str(), angle.c_str()) == 0) { m_AngleXMinimum = minimum; m_AngleXMaximum = maximum; m_AngleXMode = mode; } else if (strcmp(XML_TAG_Y_ANGLE.c_str(), angle.c_str()) == 0) { m_AngleYMinimum = minimum; m_AngleYMaximum = maximum; m_AngleYMode = mode; } } void mitk::pa::LightSource::ParseEnergy(TiXmlElement* element) { TiXmlElement* energy = element->FirstChildElement(XML_TAG_ENERGY); if (energy) { std::string energyText = energy->GetText(); m_Energy = std::stod(energyText); if (m_Verbose) std::cout << "Setting energy=" << m_Energy; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_ENERGY << "\" field in xml. Setting Energy=1" << std::endl; m_Energy = 1.0; } } void mitk::pa::LightSource::ParsePhotonSpawnArea(TiXmlElement* element) { TiXmlElement* spawnArea = element->FirstChildElement("PhotonSpawnArea"); if (spawnArea) { TiXmlElement* spawnType = spawnArea->FirstChildElement(XML_TAG_SPAWN_TYPE); if (spawnType) { std::string spawnTypeText = spawnType->GetText(); if (strcmp(XML_TAG_SPAWN_TYPE_POINT.c_str(), spawnTypeText.c_str()) == 0) { m_SpawnType = SpawnType::POINT; if (m_Verbose) std::cout << "Setting " << XML_TAG_SPAWN_TYPE << " = " << XML_TAG_SPAWN_TYPE_POINT << std::endl; } else if (strcmp(XML_TAG_SPAWN_TYPE_RECTANGLE.c_str(), spawnTypeText.c_str()) == 0) { m_SpawnType = SpawnType::RECTANGLE; if (m_Verbose) std::cout << "Setting " << XML_TAG_SPAWN_TYPE << " = " << XML_TAG_SPAWN_TYPE_RECTANGLE << std::endl; } else if (strcmp(XML_TAG_SPAWN_TYPE_CIRCLE.c_str(), spawnTypeText.c_str()) == 0) { m_SpawnType = SpawnType::CIRCLE; if (m_Verbose) std::cout << "Setting " << XML_TAG_SPAWN_TYPE << " = " << XML_TAG_SPAWN_TYPE_CIRCLE << std::endl; } else { std::cerr << "The provided SpawnType (" << spawnTypeText << ") did not match any available spawn type. Light source is not valid." << std::endl; m_IsValid = false; } } else { std::cerr << "The \"" << XML_TAG_SPAWN_TYPE << "\" element was not provided for this light source. Light source is not valid." << std::endl; m_IsValid = false; } TiXmlElement* xLocation = spawnArea->FirstChildElement(XML_TAG_X); if (xLocation) { std::string xLocationText = xLocation->GetText(); m_SpawnLocationX = std::stod(xLocationText); if (m_Verbose) std::cout << "Setting " << XML_TAG_X << "=" << m_SpawnLocationX; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_X << "\" field in xml. Setting " << XML_TAG_X << "=0" << std::endl; m_SpawnLocationX = 0; } TiXmlElement* yLocation = spawnArea->FirstChildElement(XML_TAG_Y); if (yLocation) { std::string yLocationText = yLocation->GetText(); m_SpawnLocationY = std::stod(yLocationText); if (m_Verbose) std::cout << "Setting " << XML_TAG_Y << "=" << m_SpawnLocationY; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_Y << "\" field in xml. Setting " << XML_TAG_Y << "=0" << std::endl; m_SpawnLocationY = 0; } TiXmlElement* zLocation = spawnArea->FirstChildElement(XML_TAG_Z); if (zLocation) { std::string zLocationText = zLocation->GetText(); m_SpawnLocationZ = std::stod(zLocationText); if (m_Verbose) std::cout << "Setting " << XML_TAG_Z << "=" << m_SpawnLocationZ; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_Z << "\" field in xml. Setting " << XML_TAG_Z << "=0.1" << std::endl; m_SpawnLocationZ = 0.1; } TiXmlElement* rLocation = spawnArea->FirstChildElement(XML_TAG_R); if (rLocation) { std::string rLocationText = rLocation->GetText(); m_SpawnLocationRadius = std::stod(rLocationText); if (m_Verbose) std::cout << "Setting " << XML_TAG_R << "=" << m_SpawnLocationRadius; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_R << "\" field in xml. Setting " << XML_TAG_R << "=0" << std::endl; m_SpawnLocationRadius = 0; } TiXmlElement* xLength = spawnArea->FirstChildElement(XML_TAG_X_LENGTH); if (xLength) { std::string xLengthText = xLength->GetText(); m_SpawnLocationXLength = std::stod(xLengthText); if (m_Verbose) std::cout << "Setting " << XML_TAG_X_LENGTH << "=" << m_SpawnLocationXLength << std::endl; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_X_LENGTH << "\" field in xml. Setting " << XML_TAG_X_LENGTH << "=0" << std::endl; m_SpawnLocationXLength = 0; } TiXmlElement* yLength = spawnArea->FirstChildElement(XML_TAG_Y_LENGTH); if (yLength) { std::string yLengthText = yLength->GetText(); m_SpawnLocationYLength = std::stod(yLengthText); if (m_Verbose) std::cout << "Setting " << XML_TAG_Y_LENGTH << "=" << m_SpawnLocationYLength << std::endl; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_Y_LENGTH << "\" field in xml. Setting " << XML_TAG_Y_LENGTH << "=0" << std::endl; m_SpawnLocationYLength = 0; } TiXmlElement* zLength = spawnArea->FirstChildElement(XML_TAG_Z_LENGTH); if (zLength) { std::string zLengthText = zLength->GetText(); m_SpawnLocationZLength = std::stod(zLengthText); if (m_Verbose) std::cout << "Setting " << XML_TAG_Z_LENGTH << "=" << m_SpawnLocationZLength << std::endl; } else { if (m_Verbose) std::cerr << "No \"" << XML_TAG_Z_LENGTH << "\" field in xml. Setting " << XML_TAG_Z_LENGTH << "=0" << std::endl; m_SpawnLocationZLength = 0; } } else m_IsValid = false; } mitk::pa::LightSource::PhotonInformation mitk::pa::LightSource::GetNextPhoton(double rnd1, double rnd2, double rnd3, double rnd4, double rnd5, double gau1, double gau2) { PhotonInformation returnValue; switch (m_SpawnType) { case POINT: returnValue.xPosition = m_SpawnLocationX; returnValue.yPosition = m_SpawnLocationY; returnValue.zPosition = m_SpawnLocationZ; break; case RECTANGLE: returnValue.xPosition = m_SpawnLocationX + rnd3 * m_SpawnLocationXLength; returnValue.yPosition = m_SpawnLocationY + rnd4 * m_SpawnLocationYLength; returnValue.zPosition = m_SpawnLocationZ + rnd5 * m_SpawnLocationZLength; break; case CIRCLE: double radius = rnd3 * m_SpawnLocationRadius; double angle = rnd4 * TWO_PI; returnValue.xPosition = m_SpawnLocationX + radius * cos(angle); returnValue.yPosition = m_SpawnLocationY + radius * sin(angle); returnValue.zPosition = m_SpawnLocationZ; break; } switch (m_AngleXMode) { case UNIFORM: returnValue.xAngle = rnd1 * (m_AngleXMaximum - m_AngleXMinimum) + m_AngleXMinimum; break; case GAUSSIAN: TransformResult trans = BoxMuellerTransform(gau1, gau2, (m_AngleXMaximum - m_AngleXMinimum) / 2 + m_AngleXMinimum, (m_AngleXMaximum - m_AngleXMinimum) / 2.355); returnValue.xAngle = trans.z0; break; } switch (m_AngleYMode) { case UNIFORM: returnValue.yAngle = rnd2 * (m_AngleYMaximum - m_AngleYMinimum) + m_AngleYMinimum; break; case GAUSSIAN: TransformResult trans = BoxMuellerTransform(gau1, gau2, (m_AngleYMaximum - m_AngleYMinimum) / 2 + m_AngleYMinimum, (m_AngleYMaximum - m_AngleYMinimum) / 2.355); returnValue.yAngle = trans.z1; break; } if ((returnValue.xAngle*returnValue.xAngle + returnValue.yAngle*returnValue.yAngle) > 1) { double unify = sqrt(returnValue.xAngle*returnValue.xAngle + returnValue.yAngle*returnValue.yAngle)*1.001; returnValue.xAngle = returnValue.xAngle / unify; returnValue.yAngle = returnValue.yAngle / unify; } returnValue.zAngle = sqrt(1 - returnValue.xAngle*returnValue.xAngle - returnValue.yAngle*returnValue.yAngle); if (m_Verbose) std::cout << "Created a new photon at (" << returnValue.xPosition << "|" << returnValue.yPosition << "|" << returnValue.zPosition << ") with angle (" << returnValue.xAngle << "|" << returnValue.yAngle << "|" << returnValue.zAngle << ")" << std::endl; return returnValue; } bool mitk::pa::LightSource::IsValid() { return m_IsValid; } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp index 9e9ecf4daf..03922d5835 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAPropertyCalculator.cpp @@ -1,154 +1,155 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPAPropertyCalculator.h" // us #include #include #include #include #include mitk::pa::PropertyCalculator::Properties mitk::pa::PropertyCalculator:: CalculatePropertyForSpecificWavelength( TissueType tissueType, int wavelength, double bloodOxygenInFraction) { Properties returnValue; if (!m_Valid) { mitkThrow() << "PhotoacousticPropertyGenerator was not loaded properly."; return returnValue; } double bloodOxygenation = bloodOxygenInFraction; double bloodVolumeFraction; double waterVolumeFraction; double fatVolumeFraction; double melanosomesVolumeFraction; double musp500; double fray; double bmie; double g; switch (tissueType) { case AIR: returnValue.mua = 1e-4; returnValue.mus = 1.0; returnValue.g = 1.0; return returnValue; case BLOOD: bloodVolumeFraction = 1.0; waterVolumeFraction = 0.95; fatVolumeFraction = 0.02; melanosomesVolumeFraction = 0; musp500 = 10; fray = 0; bmie = 1; g = 0.9; break; case EPIDERMIS: bloodVolumeFraction = 0; bloodOxygenation = 0; waterVolumeFraction = 0.75; fatVolumeFraction = 0.01; melanosomesVolumeFraction = 0.02; musp500 = 40; fray = 0; bmie = 1; g = 0.9; break; case FAT: bloodVolumeFraction = 0.01; bloodOxygenation = 0.75; waterVolumeFraction = 0.5; fatVolumeFraction = 0.6; melanosomesVolumeFraction = 0; musp500 = 17.47; fray = 0.2; bmie = 0.5; g = 0.9; break; case STANDARD_TISSUE: default: bloodVolumeFraction = 0.02; bloodOxygenation = 0.75; waterVolumeFraction = 0.8; fatVolumeFraction = 0.05; melanosomesVolumeFraction = 0; musp500 = 25; fray = 0.25; bmie = 1; g = 0.9; break; } // We want the reduced scattering coefficient directly. double musp = musp500 * (fray * pow(wavelength / 500.0, -4.0) + ((1 - fray) * pow(wavelength / 500.0, -bmie))); + returnValue.mus = musp; returnValue.mus = 15;//musp; double mua = bloodVolumeFraction*bloodOxygenation*m_SpectralLibMap[MapType::OXYGENATED][wavelength] + bloodVolumeFraction*(1 - bloodOxygenation)*m_SpectralLibMap[MapType::DEOXYGENATED][wavelength] + waterVolumeFraction*m_SpectralLibMap[MapType::WATER][wavelength] + fatVolumeFraction*m_SpectralLibMap[MapType::FATTY][wavelength] + melanosomesVolumeFraction*m_SpectralLibMap[MapType::MELANIN][wavelength]; returnValue.mua = mua; returnValue.g = g; return returnValue; } mitk::pa::PropertyCalculator::PropertyCalculator() { us::ModuleResource resource = us::GetModuleContext()->GetModule()->GetResource("spectralLIB.dat"); if (resource.IsValid()) { us::ModuleResourceStream resourceStream(resource); std::string line; int wavelength = 300; int counter = 0; while (std::getline(resourceStream, line, '\t')) { int currentLineIdx = counter % 6; if (currentLineIdx == 0) wavelength = std::stoi(line); else { std::istringstream lineStream(line); double tempDouble; lineStream >> tempDouble; m_SpectralLibMap[currentLineIdx][wavelength] = tempDouble; } ++counter; } } else { m_Valid = false; } m_Valid = true; } mitk::pa::PropertyCalculator::~PropertyCalculator() { m_SpectralLibMap.clear(); m_Valid = false; } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp index c2ea7abcbc..2d468a8062 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVector.cpp @@ -1,168 +1,168 @@ /*=================================================================== 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 "mitkPAVector.h" #include "chrono" #include "math.h" mitk::pa::Vector::Vector() { m_Vector.Fill(0); } mitk::pa::Vector::~Vector() { m_Vector.Fill(0); } double mitk::pa::Vector::GetNorm() { return m_Vector.GetNorm(); } double mitk::pa::Vector::GetElement(unsigned short index) { return m_Vector.GetElement(index); } void mitk::pa::Vector::SetElement(unsigned short index, double value) { m_Vector.SetElement(index, value); } void mitk::pa::Vector::Normalize() { double norm = m_Vector.GetNorm(); m_Vector.SetElement(0, m_Vector.GetElement(0) / norm); m_Vector.SetElement(1, m_Vector.GetElement(1) / norm); m_Vector.SetElement(2, m_Vector.GetElement(2) / norm); } void mitk::pa::Vector::SetValue(mitk::pa::Vector::Pointer value) { m_Vector.SetElement(0, value->GetElement(0)); m_Vector.SetElement(1, value->GetElement(1)); m_Vector.SetElement(2, value->GetElement(2)); } void mitk::pa::Vector::RandomizeByPercentage(double percentage, double bendingFactor, std::mt19937* rng) { std::uniform_real_distribution<> range(-percentage, percentage); m_Vector.SetElement(0, m_Vector.GetElement(0) + (bendingFactor * range(*rng))); m_Vector.SetElement(1, m_Vector.GetElement(1) + (bendingFactor * range(*rng))); m_Vector.SetElement(2, m_Vector.GetElement(2) + (bendingFactor * range(*rng))); } void mitk::pa::Vector::Randomize(double xLowerLimit, double xUpperLimit, double yLowerLimit, double yUpperLimit, double zLowerLimit, double zUpperLimit, std::mt19937* rng) { std::uniform_real_distribution<> rangeX(xLowerLimit, xUpperLimit); std::uniform_real_distribution<> rangeY(yLowerLimit, yUpperLimit); std::uniform_real_distribution<> rangeZ(zLowerLimit, zUpperLimit); m_Vector.SetElement(0, rangeX(*rng)); m_Vector.SetElement(1, rangeY(*rng)); m_Vector.SetElement(2, rangeZ(*rng)); } void mitk::pa::Vector::Randomize(double xLimit, double yLimit, double zLimit, std::mt19937* rng) { Randomize(0, xLimit, 0, yLimit, 0, zLimit, rng); } void mitk::pa::Vector::Randomize(std::mt19937* rng) { Randomize(-1, 1, -1, 1, -1, 1, rng); } -void mitk::pa::Vector::PrintSelf(std::ostream& os, itk::Indent indent) const +void mitk::pa::Vector::PrintSelf(std::ostream& os, itk::Indent /*indent*/) const { os << "X: " << m_Vector.GetElement(0) << std::endl; os << "Y: " << m_Vector.GetElement(1) << std::endl; os << "Z: " << m_Vector.GetElement(2) << std::endl; os << "Length: " << m_Vector.GetNorm() << std::endl; } void mitk::pa::Vector::Rotate(double thetaChange, double phiChange) { MITK_DEBUG << "Vector before rotation: (" << GetElement(0) << "|" << GetElement(1) << "|" << GetElement(2) << ")"; if (thetaChange == 0 && phiChange == 0) return; double x = GetElement(0); double y = GetElement(1); double z = GetElement(2); double r = sqrt(x*x + y*y + z*z); if (r == 0) return; double theta = acos(z / r); double phi = atan2(y, x); theta += thetaChange; phi += phiChange; SetElement(0, r * sin(theta) * cos(phi)); SetElement(1, r * sin(theta) * sin(phi)); SetElement(2, r * cos(theta)); MITK_DEBUG << "Vector after rotation: (" << GetElement(0) << "|" << GetElement(1) << "|" << GetElement(2) << ")"; } void mitk::pa::Vector::Scale(double factor) { m_Vector.SetElement(0, m_Vector.GetElement(0)*factor); m_Vector.SetElement(1, m_Vector.GetElement(1)*factor); m_Vector.SetElement(2, m_Vector.GetElement(2)*factor); } mitk::pa::Vector::Pointer mitk::pa::Vector::Clone() { auto returnVector = Vector::New(); returnVector->SetElement(0, this->GetElement(0)); returnVector->SetElement(1, this->GetElement(1)); returnVector->SetElement(2, this->GetElement(2)); return returnVector; } bool mitk::pa::Equal(const Vector::Pointer leftHandSide, const Vector::Pointer rightHandSide, double eps, bool verbose) { MITK_INFO(verbose) << "=== mitk::pa::Vector Equal ==="; if (rightHandSide.IsNull() || leftHandSide.IsNull()) { MITK_INFO(verbose) << "Cannot compare nullpointers"; return false; } if (leftHandSide->GetElement(0) - rightHandSide->GetElement(0) > eps) { MITK_INFO(verbose) << "Element[0] not equal"; return false; } if (leftHandSide->GetElement(1) - rightHandSide->GetElement(1) > eps) { MITK_INFO(verbose) << "Element[1] not equal"; return false; } if (leftHandSide->GetElement(2) - rightHandSide->GetElement(2) > eps) { MITK_INFO(verbose) << "Element[2] not equal"; return false; } return true; } diff --git a/Modules/PhotoacousticsLib/test/mitkMCThreadHandlerTest.cpp b/Modules/PhotoacousticsLib/test/mitkMCThreadHandlerTest.cpp index 675848dfdb..3a0e128fd1 100644 --- a/Modules/PhotoacousticsLib/test/mitkMCThreadHandlerTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkMCThreadHandlerTest.cpp @@ -1,115 +1,114 @@ /*=================================================================== 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 // us #include #include #include #include #include #include #include class mitkMCThreadHandlerTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkMCThreadHandlerTestSuite); MITK_TEST(testCorrectTimeMeasure); MITK_TEST(testCorrectNumberOfPhotons); MITK_TEST(testCorrectNumberOfPhotonsWithUnevenPackageSize); MITK_TEST(testCorrectNumberOfPhotonsWithTooLargePackageSize); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::MonteCarloThreadHandler::Pointer m_MonteCarloThreadHandler; long m_NumberOrTime = 500; public: void setUp() { } void testCorrectTimeMeasure() { for (int i = 0; i < 10; i++) { m_MonteCarloThreadHandler = mitk::pa::MonteCarloThreadHandler::New(m_NumberOrTime, true); auto timeBefore = std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count(); long nextWorkPackage = 0; while ((nextWorkPackage = m_MonteCarloThreadHandler->GetNextWorkPackage()) > 0) {//Do nothing } auto timeAfter = std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count(); //Assert that the time error is less than 5% in a 500ms sample size CPPUNIT_ASSERT(abs((timeAfter - timeBefore) - m_NumberOrTime) <= 25); } } void testCorrectNumberOfPhotons() { m_MonteCarloThreadHandler = mitk::pa::MonteCarloThreadHandler::New(m_NumberOrTime, false); m_MonteCarloThreadHandler->SetPackageSize(100); long numberOfPhotonsSimulated = 0; long nextWorkPackage = 0; while ((nextWorkPackage = m_MonteCarloThreadHandler->GetNextWorkPackage()) > 0) { numberOfPhotonsSimulated += nextWorkPackage; } CPPUNIT_ASSERT(numberOfPhotonsSimulated == m_NumberOrTime); } void testCorrectNumberOfPhotonsWithUnevenPackageSize() { m_MonteCarloThreadHandler = mitk::pa::MonteCarloThreadHandler::New(m_NumberOrTime, false); m_MonteCarloThreadHandler->SetPackageSize(77); long numberOfPhotonsSimulated = 0; long nextWorkPackage = 0; while ((nextWorkPackage = m_MonteCarloThreadHandler->GetNextWorkPackage()) > 0) { numberOfPhotonsSimulated += nextWorkPackage; } CPPUNIT_ASSERT(numberOfPhotonsSimulated == m_NumberOrTime); } void testCorrectNumberOfPhotonsWithTooLargePackageSize() { m_MonteCarloThreadHandler = mitk::pa::MonteCarloThreadHandler::New(m_NumberOrTime, false); m_MonteCarloThreadHandler->SetPackageSize(10000); long numberOfPhotonsSimulated = 0; long nextWorkPackage = 0; while ((nextWorkPackage = m_MonteCarloThreadHandler->GetNextWorkPackage()) > 0) { numberOfPhotonsSimulated += nextWorkPackage; } CPPUNIT_ASSERT(numberOfPhotonsSimulated == m_NumberOrTime); } void tearDown() { m_MonteCarloThreadHandler = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkMCThreadHandler) diff --git a/Modules/PhotoacousticsLib/test/mitkMcxyzXmlTest.cpp b/Modules/PhotoacousticsLib/test/mitkMcxyzXmlTest.cpp index 2c0bf44e06..7a559b0979 100644 --- a/Modules/PhotoacousticsLib/test/mitkMcxyzXmlTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkMcxyzXmlTest.cpp @@ -1,244 +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 // us #include #include #include #include #include - +#include #include #include class mitkMcxyzXmlTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkMcxyzXmlTestSuite); MITK_TEST(TestCreatePointSourceProbe); MITK_TEST(TestCreateCircleSourceProbe); MITK_TEST(TestCreateRectangleSourceProbe); MITK_TEST(TestCreateTwoPointSourcesProbe); MITK_TEST(TestCreateAllSourcesProbe); MITK_TEST(TestValuesAreInBoundsUniformRectangle); MITK_TEST(TestValuesAreInBoundsGaussianRectangle); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::LightSource::Pointer m_LightSource; mitk::pa::Probe::Pointer m_Probe; std::string m_XmlProbePointSource; std::string m_XmlProbeCircleSource; std::string m_XmlProbeRectangleSource; std::string m_XmlProbeTwoPointSources; std::string m_XmlProbeAllSources; public: void setUp() { m_LightSource = mitk::pa::LightSource::New(); LoadXmlFile("pointsource.xml", &m_XmlProbePointSource); LoadXmlFile("circlesource.xml", &m_XmlProbeCircleSource); LoadXmlFile("rectanglesource.xml", &m_XmlProbeRectangleSource); LoadXmlFile("twopointsources.xml", &m_XmlProbeTwoPointSources); LoadXmlFile("allsources.xml", &m_XmlProbeAllSources); } void LoadXmlFile(std::string filename, std::string* lines) { us::ModuleResource pointSourceXml = us::GetModuleContext()->GetModule()->GetResource(filename); std::string line; if (pointSourceXml.IsValid() && pointSourceXml.IsFile()) { us::ModuleResourceStream stream(pointSourceXml); - boost::locale::generator gen; - stream.std::istream::imbue(gen("en_GB.UTF-8")); + stream.std::istream::imbue(std::locale("C")); while (std::getline(stream, line)) { *lines = *lines + line + " "; } } else { MITK_ERROR << "Xml file was not valid"; } } void TestCreatePointSourceProbe() { m_Probe = mitk::pa::Probe::New(m_XmlProbePointSource.c_str(), true); CPPUNIT_ASSERT(true == m_Probe->IsValid()); } void TestCreateCircleSourceProbe() { m_Probe = mitk::pa::Probe::New(m_XmlProbeCircleSource.c_str(), true); CPPUNIT_ASSERT(true == m_Probe->IsValid()); } void TestCreateRectangleSourceProbe() { m_Probe = mitk::pa::Probe::New(m_XmlProbeRectangleSource.c_str(), true); CPPUNIT_ASSERT(true == m_Probe->IsValid()); } void TestCreateTwoPointSourcesProbe() { m_Probe = mitk::pa::Probe::New(m_XmlProbeTwoPointSources.c_str(), true); CPPUNIT_ASSERT(true == m_Probe->IsValid()); } void TestCreateAllSourcesProbe() { m_Probe = mitk::pa::Probe::New(m_XmlProbeAllSources.c_str(), true); CPPUNIT_ASSERT(true == m_Probe->IsValid()); } void TestValuesAreInBoundsUniformRectangle() { int MAXIMUM = 2; int MINIMUM = -2; int ANGLE_MAXIMUM = 1; int ANGLE_MINIMUM = -1; m_LightSource->SetAngleXMode(mitk::pa::LightSource::DistributionMode::UNIFORM); m_LightSource->SetAngleYMode(mitk::pa::LightSource::DistributionMode::UNIFORM); m_LightSource->SetAngleXMaximum(ANGLE_MAXIMUM); m_LightSource->SetAngleXMinimum(ANGLE_MINIMUM); m_LightSource->SetAngleYMaximum(ANGLE_MAXIMUM); m_LightSource->SetAngleYMinimum(ANGLE_MINIMUM); m_LightSource->SetSpawnLocationX(MINIMUM); m_LightSource->SetSpawnLocationXLength(2 * MAXIMUM); m_LightSource->SetSpawnLocationY(MINIMUM); m_LightSource->SetSpawnLocationYLength(2 * MAXIMUM); m_LightSource->SetSpawnLocationZ(MINIMUM); m_LightSource->SetSpawnLocationZLength(2 * MAXIMUM); m_LightSource->SetSpawnLocationRadius(MAXIMUM); m_LightSource->SetVerbose(false); m_LightSource->SetSpawnType(mitk::pa::LightSource::SpawnType::RECTANGLE); std::mt19937 rng; rng.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock:: now().time_since_epoch()).count()); std::uniform_real_distribution<> realDist(0, 1); for (int i = 0, length = 10000; i < length; i++) { auto result = m_LightSource->GetNextPhoton(realDist(rng), realDist(rng), realDist(rng), realDist(rng), realDist(rng), realDist(rng), realDist(rng)); CPPUNIT_ASSERT(result.xAngle >= ANGLE_MINIMUM); CPPUNIT_ASSERT(result.xAngle <= ANGLE_MAXIMUM); CPPUNIT_ASSERT(result.yAngle >= ANGLE_MINIMUM); CPPUNIT_ASSERT(result.yAngle <= ANGLE_MAXIMUM); CPPUNIT_ASSERT(result.zAngle >= 0); CPPUNIT_ASSERT(result.zAngle <= ANGLE_MAXIMUM); CPPUNIT_ASSERT(result.xPosition >= MINIMUM); CPPUNIT_ASSERT(result.xPosition <= MAXIMUM); CPPUNIT_ASSERT(result.yPosition >= MINIMUM); CPPUNIT_ASSERT(result.yPosition <= MAXIMUM); CPPUNIT_ASSERT(result.zPosition >= MINIMUM); CPPUNIT_ASSERT(result.zPosition <= MAXIMUM); } } void TestValuesAreInBoundsGaussianRectangle() { int MAXIMUM = 2; int MINIMUM = -2; int ANGLE_MAXIMUM = 1; int ANGLE_MINIMUM = -1; m_LightSource->SetAngleXMode(mitk::pa::LightSource::DistributionMode::GAUSSIAN); m_LightSource->SetAngleYMode(mitk::pa::LightSource::DistributionMode::GAUSSIAN); m_LightSource->SetAngleXMaximum(ANGLE_MAXIMUM); m_LightSource->SetAngleXMinimum(ANGLE_MINIMUM); m_LightSource->SetAngleYMaximum(ANGLE_MAXIMUM); m_LightSource->SetAngleYMinimum(ANGLE_MINIMUM); m_LightSource->SetSpawnLocationX(MINIMUM); m_LightSource->SetSpawnLocationXLength(2 * MAXIMUM); m_LightSource->SetSpawnLocationY(MINIMUM); m_LightSource->SetSpawnLocationYLength(2 * MAXIMUM); m_LightSource->SetSpawnLocationZ(MINIMUM); m_LightSource->SetSpawnLocationZLength(2 * MAXIMUM); m_LightSource->SetSpawnLocationRadius(MAXIMUM); m_LightSource->SetVerbose(false); m_LightSource->SetSpawnType(mitk::pa::LightSource::SpawnType::RECTANGLE); std::mt19937 rng; rng.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock:: now().time_since_epoch()).count()); std::uniform_real_distribution<> realDist(0, 1); for (int i = 0, length = 10000; i < length; i++) { auto result = m_LightSource->GetNextPhoton(realDist(rng), realDist(rng), realDist(rng), realDist(rng), realDist(rng), realDist(rng), realDist(rng)); CPPUNIT_ASSERT(result.xAngle >= ANGLE_MINIMUM); CPPUNIT_ASSERT(result.xAngle <= ANGLE_MAXIMUM); CPPUNIT_ASSERT(result.yAngle >= ANGLE_MINIMUM); CPPUNIT_ASSERT(result.yAngle <= ANGLE_MAXIMUM); CPPUNIT_ASSERT(result.zAngle >= 0); CPPUNIT_ASSERT(result.zAngle <= ANGLE_MAXIMUM); CPPUNIT_ASSERT(result.xPosition >= MINIMUM); CPPUNIT_ASSERT(result.xPosition <= MAXIMUM); CPPUNIT_ASSERT(result.yPosition >= MINIMUM); CPPUNIT_ASSERT(result.yPosition <= MAXIMUM); CPPUNIT_ASSERT(result.zPosition >= MINIMUM); CPPUNIT_ASSERT(result.zPosition <= MAXIMUM); } } void tearDown() { m_XmlProbePointSource = ""; m_XmlProbeCircleSource = ""; m_XmlProbeRectangleSource = ""; m_XmlProbeTwoPointSources = ""; m_XmlProbeAllSources = ""; m_Probe = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkMcxyzXml) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacoustic3dVolumeTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacoustic3dVolumeTest.cpp index 731ef0a328..f83ad44e06 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacoustic3dVolumeTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacoustic3dVolumeTest.cpp @@ -1,221 +1,224 @@ /*=================================================================== 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 "mitkPAVolume.h" #include #include class mitkPhotoacoustic3dVolumeTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacoustic3dVolumeTestSuite); MITK_TEST(TestCorrectGetDataAndSetDataBehavior); MITK_TEST(TestCallingConstructorWithNullParameter); MITK_TEST(TestCallingConstructorWithCorrectParameters); MITK_TEST(TestModifyImage); MITK_TEST(TestModifyComplexImage); MITK_TEST(TestConvertToMitkImage); MITK_TEST(TestDeepCopy); MITK_TEST(TestCatchException); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::Volume::Pointer m_Photoacoustic3dVolume; public: void setUp() { } void TestCallingConstructorWithNullParameter() { bool exceptionEncountered = false; try { m_Photoacoustic3dVolume = mitk::pa::Volume::New(nullptr, 3, 3, 3); } catch (...) { exceptionEncountered = true; } CPPUNIT_ASSERT(exceptionEncountered); } void TestCallingConstructorWithCorrectParameters() { double* data = new double[1]; data[0] = 3; m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetXDim() == 1); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetYDim() == 1); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetZDim() == 1); } void TestModifyImage() { double* data = new double[1]; data[0] = 3; m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); CPPUNIT_ASSERT_MESSAGE(std::to_string(m_Photoacoustic3dVolume->GetData(0, 0, 0)), m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); m_Photoacoustic3dVolume->SetData(17, 0, 0, 0); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 17); } void TestModifyComplexImage() { unsigned int xDim = 4; unsigned int yDim = 7; unsigned int zDim = 12; unsigned int length = xDim * yDim * zDim; double* data = new double[length]; for (unsigned int i = 0; i < length; i++) data[i] = 5; m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, xDim, yDim, zDim); for (unsigned int z = 0; z < zDim; z++) for (unsigned int y = 0; y < yDim; y++) for (unsigned int x = 0; x < xDim; x++) { CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(x, y, z) == 5); m_Photoacoustic3dVolume->SetData((x + y)*(z + 1), x, y, z); CPPUNIT_ASSERT(abs(m_Photoacoustic3dVolume->GetData(x, y, z) - (x + y)*(z + 1)) < mitk::eps); } } void TestCorrectGetDataAndSetDataBehavior() { unsigned int xDim = 40; unsigned int yDim = 7; unsigned int zDim = 12; unsigned int length = xDim * yDim * zDim; double* data = new double[length]; for (unsigned int i = 0; i < length; i++) data[i] = 0; m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, xDim, yDim, zDim); for (unsigned int z = 0; z < zDim; z++) for (unsigned int y = 0; y < yDim; y++) for (unsigned int x = 0; x < xDim; x++) { int index = z*xDim*yDim + x*yDim + y; m_Photoacoustic3dVolume->SetData(index, x, y, z); CPPUNIT_ASSERT_MESSAGE(std::to_string(index), m_Photoacoustic3dVolume->GetData(x, y, z) == index); } } void TestConvertToMitkImage() { double* data = new double[6]; data[0] = 3; data[1] = 3; data[2] = 3; data[3] = 3; data[4] = 3; data[5] = 3; m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 2, 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 2) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 0) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 2) == 3); m_Photoacoustic3dVolume->SetData(17, 0, 0, 0); m_Photoacoustic3dVolume->SetData(17, 0, 1, 0); m_Photoacoustic3dVolume->SetData(17, 0, 1, 2); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 17); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 2) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 0) == 17); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 1) == 3); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 1, 2) == 17); mitk::Image::Pointer mitkImage = m_Photoacoustic3dVolume->AsMitkImage(); CPPUNIT_ASSERT(mitkImage->GetDimensions()[0] == 2); CPPUNIT_ASSERT(mitkImage->GetDimensions()[1] == 1); CPPUNIT_ASSERT(mitkImage->GetDimensions()[2] == 3); mitk::ImageReadAccessor readAccess(mitkImage, mitkImage->GetVolumeData()); double* copyData = (double*)readAccess.GetData(); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[0]), copyData[0] == 17); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[1]), copyData[1] == 17); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[2]), copyData[2] == 3); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[3]), copyData[3] == 3); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[4]), copyData[4] == 3); CPPUNIT_ASSERT_MESSAGE(std::to_string(copyData[5]), copyData[5] == 17); } void TestDeepCopy() { double* data = new double[1]; data[0] = 3; m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); mitk::pa::Volume::Pointer copiedVolume = m_Photoacoustic3dVolume->DeepCopy(); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetXDim() == copiedVolume->GetXDim()); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetYDim() == copiedVolume->GetYDim()); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetZDim() == copiedVolume->GetZDim()); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 3); CPPUNIT_ASSERT(copiedVolume->GetData(0, 0, 0) == 3); m_Photoacoustic3dVolume->SetData(17, 0, 0, 0); CPPUNIT_ASSERT(m_Photoacoustic3dVolume->GetData(0, 0, 0) == 17); CPPUNIT_ASSERT(copiedVolume->GetData(0, 0, 0) == 3); } void AssertIndexException(unsigned int x, unsigned int y, unsigned int z) { bool exceptionCaught = false; try { double thisIsIrrelevant = m_Photoacoustic3dVolume->GetData(x, y, z); + thisIsIrrelevant += 1; } catch (...) { exceptionCaught = true; + if (exceptionCaught) + exceptionCaught = true; } #ifdef _DEBUG CPPUNIT_ASSERT(exceptionCaught); #endif - } + } void TestCatchException() { double* data = new double[1]; data[0] = 3; m_Photoacoustic3dVolume = mitk::pa::Volume::New(data, 1, 1, 1); AssertIndexException(1, 0, 0); AssertIndexException(0, 1, 0); AssertIndexException(0, 0, 1); AssertIndexException(18, 1, 222); } void tearDown() { m_Photoacoustic3dVolume = nullptr; } -}; + }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacoustic3dVolume) diff --git a/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp index f44cc61851..37d662697f 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticIOTest.cpp @@ -1,183 +1,187 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include #include #include #include #include #include #include #include -#include +#include 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); m_Test3DVolume = createTest3DVolume(5); - - CPPUNIT_ASSERT(boost::filesystem::create_directories(TEST_FOLDER_PATH)); - CPPUNIT_ASSERT(boost::filesystem::create_directories(TEST_QUALIFIED_FOLDER_PATH)); - CPPUNIT_ASSERT(boost::filesystem::create_directories(TEST_FOLDER_PATH + FOLDER_FOLDER + FOLDER_FOLDER)); - CPPUNIT_ASSERT(boost::filesystem::create_directories(FCM_PATH)); + 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); } mitk::pa::TissueGeneratorParameters::Pointer createTestVolumeParameters() { auto returnParameters = mitk::pa::TissueGeneratorParameters::New(); returnParameters->SetXDim(10); returnParameters->SetYDim(10); returnParameters->SetZDim(10); returnParameters->SetBackgroundAbsorption(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..", boost::filesystem::remove_all(TEST_FOLDER_PATH) > 0); + //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/mitkPhotoacousticVectorTest.cpp b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVectorTest.cpp index ea90456a56..4cbe2b9f1b 100644 --- a/Modules/PhotoacousticsLib/test/mitkPhotoacousticVectorTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkPhotoacousticVectorTest.cpp @@ -1,267 +1,260 @@ /*=================================================================== 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 #define _USE_MATH_DEFINES #include #include "mitkPAVector.h" class mitkPhotoacousticVectorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkPhotoacousticVectorTestSuite); MITK_TEST(TestNormalizeVector); MITK_TEST(TestRotateVectorZeroDegrees); MITK_TEST(TestRotatedVectorPositiveDegrees); MITK_TEST(TestRotateVectorZeroDegrees); MITK_TEST(TestScaleVector); MITK_TEST(TestCloneVector); CPPUNIT_TEST_SUITE_END(); private: mitk::pa::Vector::Pointer m_TestVector; mitk::pa::Vector::Pointer m_TestReturnVector; const double DIF_VAL = 0.001; const double TWO_PI = 6.283185; public: void setUp() { m_TestVector = mitk::pa::Vector::New(); m_TestReturnVector = mitk::pa::Vector::New(); } void TestNormalizeVector() { std::stringstream output; int a = 2; int b = 3; int c = 4; m_TestVector->SetElement(0, a); m_TestVector->SetElement(1, b); m_TestVector->SetElement(2, c); output << "The vectorlength should be"; output << sqrt(a*a + b*b + c*c); CPPUNIT_ASSERT_EQUAL_MESSAGE(output.str(), sqrt(a*a + b*b + c*c), m_TestVector->GetNorm()); output.flush(); m_TestVector->Normalize(); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vectorlength should be 1.", true, m_TestVector->GetNorm() - 1 < DIF_VAL); } void TestRotateVectorZeroDegrees() { int a = 1; int b = 2; int c = 3; double length; m_TestVector->SetElement(0, a); m_TestVector->SetElement(1, b); m_TestVector->SetElement(2, c); length = m_TestVector->GetNorm(); m_TestVector->Rotate(0, 0); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", length, m_TestVector->GetNorm()); CPPUNIT_ASSERT_MESSAGE("The vector value at index0 should be 1.0", m_TestVector->GetElement(0) - 1 < DIF_VAL); CPPUNIT_ASSERT_MESSAGE("The vector value at index1 should be 2.0", m_TestVector->GetElement(1) - 2 < DIF_VAL); CPPUNIT_ASSERT_MESSAGE("The vector value at index2 should be 3.0", m_TestVector->GetElement(2) - 3 < DIF_VAL); } void TestRotatedVectorPositiveDegrees() { MITK_INFO << atan2(0, 0); for (int r = 0; r < 10; r++) { for (double phi = 0.1; phi < 3; phi += 0.1) { for (double theta = 0.1; theta < 3; theta += 0.1) { double rotateTheta = 0.1; double rotatePhi = 0.1; m_TestVector->SetElement(0, r * sin(theta) * cos(phi)); m_TestVector->SetElement(1, r * sin(theta) * sin(phi)); m_TestVector->SetElement(2, r * cos(theta)); m_TestVector->Rotate(rotateTheta, rotatePhi); double newTheta = fmod(theta + rotateTheta, TWO_PI); double newPhi = fmod(phi + rotatePhi, TWO_PI); double expectedX = r * sin(newTheta) * cos(newPhi); double expectedY = r * sin(newTheta) * sin(newPhi); double expectedZ = r * cos(newTheta); CPPUNIT_ASSERT_MESSAGE("The vector value at index0 should be " + std::to_string(expectedX) + " but was " + std::to_string(m_TestVector->GetElement(0)) + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta), m_TestVector->GetElement(0) - expectedX < DIF_VAL); CPPUNIT_ASSERT_MESSAGE("The vector value at index1 should be " + std::to_string(expectedY) + " but was " + std::to_string(m_TestVector->GetElement(0)) + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta), m_TestVector->GetElement(1) - expectedY < DIF_VAL); CPPUNIT_ASSERT_MESSAGE("The vector value at index2 should be " + std::to_string(expectedZ) + " but was " + std::to_string(m_TestVector->GetElement(0)) + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta), m_TestVector->GetElement(2) - expectedZ < DIF_VAL); } } } } void TestRotatedVectorNegativeDegrees() { for (int r = 0; r < 10; r++) { for (double phi = -0.1; phi > -3; phi -= 0.1) { for (double theta = -0.1; theta > -3; theta -= 0.1) { double rotateTheta = -0.1; double rotatePhi = -0.1; m_TestVector->SetElement(0, r * sin(theta) * cos(phi)); m_TestVector->SetElement(1, r * sin(theta) * sin(phi)); m_TestVector->SetElement(2, r * cos(theta)); m_TestVector->Rotate(rotateTheta, rotatePhi); double newTheta = fmod(theta + rotateTheta, TWO_PI); double newPhi = fmod(phi + rotatePhi, TWO_PI); double expectedX = r * sin(newTheta) * cos(newPhi); double expectedY = r * sin(newTheta) * sin(newPhi); double expectedZ = r * cos(newTheta); CPPUNIT_ASSERT_MESSAGE("The vector value at index0 should be " + std::to_string(expectedX) + " but was " + std::to_string(m_TestVector->GetElement(0)) + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta), m_TestVector->GetElement(0) - expectedX < DIF_VAL); CPPUNIT_ASSERT_MESSAGE("The vector value at index1 should be " + std::to_string(expectedY) + " but was " + std::to_string(m_TestVector->GetElement(0)) + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta), m_TestVector->GetElement(1) - expectedY < DIF_VAL); CPPUNIT_ASSERT_MESSAGE("The vector value at index2 should be " + std::to_string(expectedZ) + " but was " + std::to_string(m_TestVector->GetElement(0)) + " at r=" + std::to_string(r) + " phi=" + std::to_string(phi) + " theta=" + std::to_string(theta), m_TestVector->GetElement(2) - expectedZ < DIF_VAL); } } } } void TestScaleVector() { double a = 1.0; double b = 2.0; double c = 3.0; - double length; - for (double testFactor = -2.0; testFactor <= 2.0; testFactor += 0.3) { double potElement0Fctr; double potElement1Fctr; double potElement2Fctr; - double testLength; std::stringstream output; m_TestVector->SetElement(0, a); m_TestVector->SetElement(1, b); m_TestVector->SetElement(2, c); - length = m_TestVector->GetNorm(); - potElement0Fctr = (m_TestVector->GetElement(0)*testFactor)*(m_TestVector->GetElement(0)*testFactor); potElement1Fctr = (m_TestVector->GetElement(1)*testFactor)*(m_TestVector->GetElement(1)*testFactor); potElement2Fctr = (m_TestVector->GetElement(2)*testFactor)*(m_TestVector->GetElement(2)*testFactor); - testLength = sqrt(potElement0Fctr + potElement1Fctr + potElement2Fctr); - m_TestVector->Scale(testFactor); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should not be equal", sqrt(potElement0Fctr + potElement1Fctr + potElement2Fctr), m_TestVector->GetNorm()); output << "The vector value at index0 should be"; output << a*testFactor; CPPUNIT_ASSERT_EQUAL_MESSAGE(output.str(), a*testFactor, m_TestVector->GetElement(0)); output.flush(); output << "The vector value at index1 should be"; output << b*testFactor; CPPUNIT_ASSERT_EQUAL_MESSAGE(output.str(), b*testFactor, m_TestVector->GetElement(1)); output.flush(); output << "The vector value at index2 should be"; output << c*testFactor; CPPUNIT_ASSERT_EQUAL_MESSAGE(output.str(), c*testFactor, m_TestVector->GetElement(2)); output.flush(); } } void TestCloneVector() { int a = 1; int b = 2; int c = 3; m_TestVector->SetElement(0, a); m_TestVector->SetElement(1, b); m_TestVector->SetElement(2, c); m_TestReturnVector = m_TestVector->Clone(); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector length should be equal", (m_TestVector->GetNorm()), m_TestReturnVector->GetNorm()); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be equal", m_TestVector->GetElement(0), m_TestReturnVector->GetElement(0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index1 should be equal", m_TestVector->GetElement(1), m_TestReturnVector->GetElement(1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index2 should be equal", m_TestVector->GetElement(2), m_TestReturnVector->GetElement(2)); m_TestReturnVector->Rotate(M_PI / 4, M_PI / 4); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be not equal", true, m_TestVector->GetElement(0) != m_TestReturnVector->GetElement(0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be not equal", true, m_TestVector->GetElement(1) != m_TestReturnVector->GetElement(1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be not equal", true, m_TestVector->GetElement(2) != m_TestReturnVector->GetElement(2)); for (double testFactor = -2.0; testFactor <= 2.0; testFactor += 0.3) { m_TestReturnVector->Scale(testFactor); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be not equal", true, m_TestVector->GetElement(0) != m_TestReturnVector->GetElement(0)); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be not equal", true, m_TestVector->GetElement(1) != m_TestReturnVector->GetElement(1)); CPPUNIT_ASSERT_EQUAL_MESSAGE("The vector value at index0 should be not equal", true, m_TestVector->GetElement(2) != m_TestReturnVector->GetElement(2)); } } void tearDown() { m_TestVector = nullptr; m_TestReturnVector = nullptr; } }; MITK_TEST_SUITE_REGISTRATION(mitkPhotoacousticVector) diff --git a/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp b/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp index 4d8f1613e4..305057568e 100644 --- a/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp +++ b/Modules/PhotoacousticsLib/test/mitkSimulationBatchGeneratorTest.cpp @@ -1,89 +1,90 @@ ///*=================================================================== //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 class mitkSimulationBatchGeneratorTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkSimulationBatchGeneratorTestSuite); MITK_TEST(testGenerateBatchFileString); MITK_TEST(testGenerateBatchFileAndSaveFile); CPPUNIT_TEST_SUITE_END(); private: const std::string TEST_FOLDER_PATH = "testFiles/"; mitk::pa::SimulationBatchGeneratorParameters::Pointer m_Parameters; mitk::pa::Volume::Pointer m_Test3DVolume; public: void setUp() { m_Parameters = mitk::pa::SimulationBatchGeneratorParameters::New(); m_Parameters->SetBinaryPath("binary"); m_Parameters->SetNrrdFilePath(TEST_FOLDER_PATH); m_Parameters->SetNumberOfPhotons(100); m_Parameters->SetTissueName("tissueName"); m_Parameters->SetVolumeIndex(0); m_Parameters->SetYOffsetLowerThresholdInCentimeters(-1); m_Parameters->SetYOffsetUpperThresholdInCentimeters(1); m_Parameters->SetYOffsetStepInCentimeters(0.5); m_Test3DVolume = createTest3DVolume(5); - CPPUNIT_ASSERT(boost::filesystem::create_directories(TEST_FOLDER_PATH)); + itk::FileTools::CreateDirectory(TEST_FOLDER_PATH); + CPPUNIT_ASSERT(itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH)); } 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); } void testGenerateBatchFileString() { std::string batchGenerationString = mitk::pa::SimulationBatchGenerator::CreateBatchSimulationString(m_Parameters); CPPUNIT_ASSERT(!batchGenerationString.empty()); } void testGenerateBatchFileAndSaveFile() { mitk::pa::SimulationBatchGenerator::WriteBatchFileAndSaveTissueVolume(m_Parameters, m_Test3DVolume->AsMitkImage()); - CPPUNIT_ASSERT(boost::filesystem::exists(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000.nrrd")); - CPPUNIT_ASSERT(boost::filesystem::exists(TEST_FOLDER_PATH + "simulate_all.sh") || boost::filesystem::exists(TEST_FOLDER_PATH + "simulate_all.bat")); - CPPUNIT_ASSERT(boost::filesystem::exists(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000") - && boost::filesystem::is_directory(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000")); + CPPUNIT_ASSERT(itksys::SystemTools::FileExists(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000.nrrd")); + CPPUNIT_ASSERT(itksys::SystemTools::FileExists(TEST_FOLDER_PATH + "simulate_all.sh") || itksys::SystemTools::FileExists(TEST_FOLDER_PATH + "simulate_all.bat")); + CPPUNIT_ASSERT(itksys::SystemTools::FileExists(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000") + && itksys::SystemTools::FileIsDirectory(TEST_FOLDER_PATH + m_Parameters->GetTissueName() + "000")); } void tearDown() { m_Parameters = nullptr; - CPPUNIT_ASSERT_MESSAGE("Resource leak of test files onto hard drive..", boost::filesystem::remove_all(TEST_FOLDER_PATH) > 0); + CPPUNIT_ASSERT_MESSAGE("Resource leak of test files onto hard drive..", itksys::SystemTools::RemoveADirectory(TEST_FOLDER_PATH) == true); } }; MITK_TEST_SUITE_REGISTRATION(mitkSimulationBatchGenerator) diff --git a/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulator.cpp b/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulator.cpp index f582116d87..e1b432613a 100644 --- a/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulator.cpp +++ b/Plugins/org.mitk.gui.qt.photoacoustics.simulation/src/internal/PASimulator.cpp @@ -1,271 +1,271 @@ /*=================================================================== 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. ===================================================================*/ // Blueberry #include #include // Qmitk #include "PASimulator.h" // Qt #include #include #include // mitk #include #include #include #include #include const std::string PASimulator::VIEW_ID = "org.mitk.views.pasimulator"; void PASimulator::SetFocus() { m_Controls.pushButtonShowRandomTissue->setFocus(); } void PASimulator::CreateQtPartControl(QWidget *parent) { m_Controls.setupUi(parent); connect(m_Controls.pushButtonShowRandomTissue, SIGNAL(clicked()), this, SLOT(DoImageProcessing())); connect(m_Controls.checkBoxGauss, SIGNAL(stateChanged(int)), this, SLOT(ClickedGaussBox())); connect(m_Controls.pushButtonOpenPath, SIGNAL(clicked()), this, SLOT(OpenFolder())); connect(m_Controls.pushButtonOpenBinary, SIGNAL(clicked()), this, SLOT(OpenBinary())); connect(m_Controls.checkBoxGenerateBatch, SIGNAL(clicked()), this, SLOT(UpdateVisibilityOfBatchCreation())); connect(m_Controls.pushButtonAjustWavelength, SIGNAL(clicked()), this, SLOT(UpdateParametersAccordingToWavelength())); connect(m_Controls.checkBoxRngSeed, SIGNAL(clicked()), this, SLOT(ClickedCheckboxFixedSeed())); connect(m_Controls.checkBoxRandomizeParameters, SIGNAL(clicked()), this, SLOT(ClickedRandomizePhysicalParameters())); m_Controls.spinboxSigma->setEnabled(false); m_Controls.labelSigma->setEnabled(false); - char* home_env = std::getenv("HOME"); - if (home_env == nullptr) + std::string home_env = std::string(std::getenv("HOME")); + if (home_env.empty()) { - home_env = std::getenv("HOMEPATH"); + home_env = std::string(std::getenv("HOMEPATH")); } - if (home_env == nullptr) + if (home_env.empty()) { home_env = ""; } - m_Controls.label_NrrdFilePath->setText(home_env); + m_Controls.label_NrrdFilePath->setText(home_env.c_str()); m_PhotoacousticPropertyCalculator = mitk::pa::PropertyCalculator::New(); UpdateVisibilityOfBatchCreation(); ClickedRandomizePhysicalParameters(); ClickedCheckboxFixedSeed(); ClickedGaussBox(); } void PASimulator::ClickedRandomizePhysicalParameters() { m_Controls.spinboxRandomizeParameters->setEnabled(m_Controls.checkBoxRandomizeParameters->isChecked()); } void PASimulator::ClickedCheckboxFixedSeed() { m_Controls.spinBoxRngSeed->setEnabled(m_Controls.checkBoxRngSeed->isChecked()); } void PASimulator::UpdateParametersAccordingToWavelength() { int wavelength = m_Controls.spinboxWavelength->value(); double bloodOxygenation = m_Controls.spinboxBloodOxygenSaturation->value() / 100; auto result = m_PhotoacousticPropertyCalculator->CalculatePropertyForSpecificWavelength( mitk::pa::PropertyCalculator::TissueType::BLOOD, wavelength, bloodOxygenation); m_Controls.spinboxMaxAbsorption->setValue(result.mua); m_Controls.spinboxMinAbsorption->setValue(result.mua); m_Controls.spinboxBloodVesselScatteringMinimum->setValue(result.mus); m_Controls.spinboxBloodVesselScatteringMaximum->setValue(result.mus); m_Controls.spinboxBloodVesselAnisotropyMinimum->setValue(result.g); m_Controls.spinboxBloodVesselAnisotropyMaximum->setValue(result.g); result = m_PhotoacousticPropertyCalculator->CalculatePropertyForSpecificWavelength( mitk::pa::PropertyCalculator::TissueType::EPIDERMIS, wavelength, bloodOxygenation); m_Controls.spinboxSkinAbsorption->setValue(result.mua); m_Controls.spinboxSkinScattering->setValue(result.mus); m_Controls.spinboxSkinAnisotropy->setValue(result.g); result = m_PhotoacousticPropertyCalculator->CalculatePropertyForSpecificWavelength( mitk::pa::PropertyCalculator::TissueType::STANDARD_TISSUE, wavelength, bloodOxygenation); m_Controls.spinboxBackgroundAbsorption->setValue(result.mua); m_Controls.spinboxBackgroundScattering->setValue(result.mus); m_Controls.spinboxBackgroundAnisotropy->setValue(result.g); } void PASimulator::UpdateVisibilityOfBatchCreation() { m_Controls.widgetBatchFile->setVisible(m_Controls.checkBoxGenerateBatch->isChecked()); } mitk::pa::TissueGeneratorParameters::Pointer PASimulator::GetParametersFromUIInput() { auto parameters = mitk::pa::TissueGeneratorParameters::New(); // Getting settings from UI // General settings parameters->SetXDim(m_Controls.spinboxXDim->value()); parameters->SetYDim(m_Controls.spinboxYDim->value()); parameters->SetZDim(m_Controls.spinboxZDim->value()); parameters->SetDoVolumeSmoothing(m_Controls.checkBoxGauss->isChecked()); if (parameters->GetDoVolumeSmoothing()) parameters->SetVolumeSmoothingSigma(m_Controls.spinboxSigma->value()); parameters->SetRandomizePhysicalProperties(m_Controls.checkBoxRandomizeParameters->isChecked()); parameters->SetRandomizePhysicalPropertiesPercentage(m_Controls.spinboxRandomizeParameters->value()); parameters->SetVoxelSpacingInCentimeters(m_Controls.spinboxSpacing->value()); parameters->SetUseRngSeed(m_Controls.checkBoxRngSeed->isChecked()); parameters->SetRngSeed(m_Controls.spinBoxRngSeed->value()); // Monte Carlo simulation parameters parameters->SetMCflag(m_Controls.spinboxMcFlag->value()); parameters->SetMCLaunchflag(m_Controls.spinboxLaunchFlag->value()); parameters->SetMCBoundaryflag(m_Controls.spinboxboundaryFlag->value()); parameters->SetMCLaunchPointX(m_Controls.spinboxLaunchpointX->value()); parameters->SetMCLaunchPointY(m_Controls.spinboxLaunchpointY->value()); parameters->SetMCLaunchPointZ(m_Controls.spinboxLaunchpointZ->value()); parameters->SetMCFocusPointX(m_Controls.spinboxFocuspointX->value()); parameters->SetMCFocusPointY(m_Controls.spinboxFocuspointY->value()); parameters->SetMCFocusPointZ(m_Controls.spinboxFocuspointZ->value()); parameters->SetMCTrajectoryVectorX(m_Controls.spinboxTrajectoryVectorX->value()); parameters->SetMCTrajectoryVectorY(m_Controls.spinboxTrajectoryVectorY->value()); parameters->SetMCTrajectoryVectorZ(m_Controls.spinboxTrajectoryVectorZ->value()); parameters->SetMCRadius(m_Controls.spinboxRadius->value()); parameters->SetMCWaist(m_Controls.spinboxWaist->value()); // Vessel settings parameters->SetMaxVesselAbsorption(m_Controls.spinboxMaxAbsorption->value()); parameters->SetMinVesselAbsorption(m_Controls.spinboxMinAbsorption->value()); parameters->SetMaxVesselBending(m_Controls.spinboxMaxBending->value()); parameters->SetMinVesselBending(m_Controls.spinboxMinBending->value()); parameters->SetMaxVesselRadiusInMillimeters(m_Controls.spinboxMaxDiameter->value()); parameters->SetMinVesselRadiusInMillimeters(m_Controls.spinboxMinDiameter->value()); parameters->SetMaxNumberOfVessels(m_Controls.spinboxMaxVessels->value()); parameters->SetMinNumberOfVessels(m_Controls.spinboxMinVessels->value()); parameters->SetMinVesselScattering(m_Controls.spinboxBloodVesselScatteringMinimum->value()); parameters->SetMaxVesselScattering(m_Controls.spinboxBloodVesselScatteringMaximum->value()); parameters->SetMinVesselAnisotropy(m_Controls.spinboxBloodVesselAnisotropyMinimum->value()); parameters->SetMaxVesselAnisotropy(m_Controls.spinboxBloodVesselAnisotropyMaximum->value()); parameters->SetVesselBifurcationFrequency(m_Controls.spinboxBifurcationFrequency->value()); parameters->SetMinVesselZOrigin(m_Controls.spinboxMinSpawnDepth->value()); parameters->SetMaxVesselZOrigin(m_Controls.spinboxMaxSpawnDepth->value()); // Background tissue settings parameters->SetBackgroundAbsorption(m_Controls.spinboxBackgroundAbsorption->value()); parameters->SetBackgroundScattering(m_Controls.spinboxBackgroundScattering->value()); parameters->SetBackgroundAnisotropy(m_Controls.spinboxBackgroundAnisotropy->value()); // Air settings parameters->SetAirThicknessInMillimeters(m_Controls.spinboxAirThickness->value()); //Skin tissue settings parameters->SetSkinThicknessInMillimeters(m_Controls.spinboxSkinThickness->value()); parameters->SetSkinAbsorption(m_Controls.spinboxSkinAbsorption->value()); parameters->SetSkinScattering(m_Controls.spinboxSkinScattering->value()); parameters->SetSkinAnisotropy(m_Controls.spinboxSkinAnisotropy->value()); parameters->SetCalculateNewVesselPositionCallback(&mitk::pa::VesselMeanderStrategy::CalculateRandomlyDivergingPosition); return parameters; } void PASimulator::DoImageProcessing() { int numberOfVolumes = 1; if (m_Controls.checkBoxGenerateBatch->isChecked()) { if (m_Controls.labelBinarypath->text().isNull() || m_Controls.labelBinarypath->text().isEmpty()) { QMessageBox::warning(nullptr, QString("Warning"), QString("You need to specify the binary first!")); return; } numberOfVolumes = m_Controls.spinboxNumberVolumes->value(); if (numberOfVolumes < 1) { QMessageBox::warning(nullptr, QString("Warning"), QString("You need to create at least one volume!")); return; } } auto tissueParameters = GetParametersFromUIInput(); for (int volumeIndex = 0; volumeIndex < numberOfVolumes; volumeIndex++) { mitk::pa::InSilicoTissueVolume::Pointer volume = mitk::pa::InSilicoTissueGenerator::GenerateInSilicoData(tissueParameters); mitk::Image::Pointer tissueVolume = volume->ConvertToMitkImage(); if (m_Controls.checkBoxGenerateBatch->isChecked()) { std::string nrrdFilePath = m_Controls.label_NrrdFilePath->text().toStdString(); std::string tissueName = m_Controls.lineEditTissueName->text().toStdString(); std::string binaryPath = m_Controls.labelBinarypath->text().toStdString(); long numberOfPhotons = m_Controls.spinboxNumberPhotons->value() * 1000L; auto batchParameters = mitk::pa::SimulationBatchGeneratorParameters::New(); batchParameters->SetBinaryPath(binaryPath); batchParameters->SetNrrdFilePath(nrrdFilePath); batchParameters->SetNumberOfPhotons(numberOfPhotons); batchParameters->SetTissueName(tissueName); batchParameters->SetVolumeIndex(volumeIndex); batchParameters->SetYOffsetLowerThresholdInCentimeters(m_Controls.spinboxFromValue->value()); batchParameters->SetYOffsetUpperThresholdInCentimeters(m_Controls.spinboxToValue->value()); batchParameters->SetYOffsetStepInCentimeters(m_Controls.spinboxStepValue->value()); mitk::pa::SimulationBatchGenerator::WriteBatchFileAndSaveTissueVolume(batchParameters, tissueVolume); } else { mitk::DataNode::Pointer dataNode = mitk::DataNode::New(); dataNode->SetData(tissueVolume); dataNode->SetName(m_Controls.lineEditTissueName->text().toStdString()); this->GetDataStorage()->Add(dataNode); mitk::RenderingManager::GetInstance()->InitializeViewsByBoundingObjects(this->GetDataStorage()); } } } void PASimulator::ClickedGaussBox() { if (m_Controls.checkBoxGauss->isChecked()) { m_Controls.spinboxSigma->setEnabled(true); m_Controls.labelSigma->setEnabled(true); } else { m_Controls.spinboxSigma->setEnabled(false); m_Controls.labelSigma->setEnabled(false); } } void PASimulator::OpenFolder() { m_Controls.label_NrrdFilePath->setText(QFileDialog::getExistingDirectory().append("/")); } void PASimulator::OpenBinary() { m_Controls.labelBinarypath->setText(QFileDialog::getOpenFileName()); }