diff --git a/Modules/PhotoacousticsLib/CMakeLists.txt b/Modules/PhotoacousticsLib/CMakeLists.txt index c6361e065a..8c89b2658e 100644 --- a/Modules/PhotoacousticsLib/CMakeLists.txt +++ b/Modules/PhotoacousticsLib/CMakeLists.txt @@ -1,12 +1,12 @@ MITK_CREATE_MODULE( INCLUDE_DIRS PUBLIC include INTERNAL_INCLUDE_DIRS ${INCLUDE_DIRS_INTERNAL} DEPENDS PUBLIC MitkAlgorithmsExt tinyxml PACKAGE_DEPENDS tinyxml PUBLIC ITK ) -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..c6ab1d8e5f 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Volume/mitkPAInSilicoTissueVolume.cpp @@ -1,317 +1,321 @@ /*=================================================================== 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/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/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/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)