diff --git a/Modules/PhotoacousticSimulation/CMakeLists.txt b/Modules/PhotoacousticSimulation/CMakeLists.txt index cad1acd1c2..096d199cd1 100644 --- a/Modules/PhotoacousticSimulation/CMakeLists.txt +++ b/Modules/PhotoacousticSimulation/CMakeLists.txt @@ -1,10 +1,11 @@ MITK_CREATE_MODULE( PhotoacousticSimulation - INCLUDE_DIRS Algorithms Utils - DEPENDS MitkAlgorithmsExt + INCLUDE_DIRS src/Algorithms src/Utils src/ProbeDesign + DEPENDS MitkAlgorithmsExt tinyxml PACKAGE_DEPENDS + tinyxml PUBLIC ITK ) add_subdirectory(MCxyz) add_subdirectory(Testing) diff --git a/Modules/PhotoacousticSimulation/MCxyz/CMakeLists.txt b/Modules/PhotoacousticSimulation/MCxyz/CMakeLists.txt index 4e90d3b80c..c5b601ab7b 100644 --- a/Modules/PhotoacousticSimulation/MCxyz/CMakeLists.txt +++ b/Modules/PhotoacousticSimulation/MCxyz/CMakeLists.txt @@ -1,9 +1,9 @@ OPTION(BUILD_PhotoacousticSimulationMCxyz "Build MiniApp for Photoacoustic Simulation Module (mcxyz)" ON) IF(BUILD_PhotoacousticSimulationMCxyz) PROJECT( MCxyz ) mitk_create_executable(MCxyz - DEPENDS MitkCommandLine MitkCore #mitkLog mitkLogMacros #MitkPhotoacousticSimulation + DEPENDS MitkCommandLine MitkCore MitkPhotoacousticSimulation PACKAGE_DEPENDS CPP_FILES mitkMCxyz.cpp) ENDIF() diff --git a/Modules/PhotoacousticSimulation/MCxyz/files.cmake b/Modules/PhotoacousticSimulation/MCxyz/files.cmake index ca6f6a067c..91f4bbcb81 100644 --- a/Modules/PhotoacousticSimulation/MCxyz/files.cmake +++ b/Modules/PhotoacousticSimulation/MCxyz/files.cmake @@ -1,3 +1,3 @@ set(CPP_FILES - mitkMCxyz.cpp + mitkMCxyz.cpp ) diff --git a/Modules/PhotoacousticSimulation/MCxyz/mitkMCxyz.cpp b/Modules/PhotoacousticSimulation/MCxyz/mitkMCxyz.cpp index dc21649633..a8dd7b1594 100755 --- a/Modules/PhotoacousticSimulation/MCxyz/mitkMCxyz.cpp +++ b/Modules/PhotoacousticSimulation/MCxyz/mitkMCxyz.cpp @@ -1,1515 +1,1565 @@ /*=================================================================== 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 "mitkPhotoacousticLightSource.h" + #ifdef __linux__ #include #include #else #include #endif #define Ntiss 19 /* Number of tissue types. */ #define ls 1.0E-7 /* Moving photon a little bit off the voxel face */ #define PI 3.1415926 #define LIGHTSPEED 2.997925E10 /* in vacuo speed of light [cm/s] */ #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 COS90D 1.0E-6 /* If cos(theta) <= COS90D, theta >= PI/2 - 1e-6 rad. */ #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; DetectorVoxel(Location location, long totalNumberOfVoxels) { this->location = location; this->fluenceContribution = (double *)malloc(totalNumberOfVoxels*sizeof(double)); for(int j=0; j* DetectorVoxelList = new std::vector(); mitk::Image::Pointer m_inputImage; InputValues() { tissueType = nullptr; muaVector = nullptr; musVector = nullptr; gVector = nullptr; mcflag = 0; launchflag = 0; boundaryflag = 0; } void LoadValues (std::string localInputFilename, float yOffset) { inputFilename = localInputFilename; try { m_inputImage = mitk::IOUtil::LoadImage(inputFilename); } catch(...) { if(verbose) std::cout << "No .nrrd file found ... switching to legacy mode." << std::endl; } 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(); return; } /* IMPORT myname_H.mci */ headerFilename = inputFilename + "_H.mci"; fid = fopen(headerFilename.c_str(),"r"); if(fid == nullptr) { std::cout << "Could not open file " << headerFilename << ". Please check the arguments." << std::endl; exit(EXIT_FAILURE); } fgets(buf, 32, fid); // run parameters sscanf(buf, "%lf", &simulationTimeFromFile); // desired time duration of run [min] fgets(buf, 32, fid); sscanf(buf, "%d", &Nx); // # of bins fgets(buf, 32,fid); sscanf(buf, "%d", &Ny); // # of bins fgets(buf, 32,fid); sscanf(buf, "%d", &Nz); // # of bins fgets(buf, 32,fid); sscanf(buf, "%lf", &xSpacing); // size of bins [cm] fgets(buf, 32,fid); sscanf(buf, "%lf", &ySpacing); // size of bins [cm] fgets(buf, 32,fid); sscanf(buf, "%lf", &zSpacing); // size of bins [cm] // launch parameters fgets(buf, 32,fid); sscanf(buf, "%d", &mcflag); // mcflag, 0 = uniform, 1 = Gaussian, 2 = iso-pt fgets(buf, 32,fid); sscanf(buf, "%d", &launchflag); // launchflag, 0 = ignore, 1 = manually set fgets(buf, 32,fid); sscanf(buf, "%d", &boundaryflag); // 0 = no boundaries, 1 = escape at all boundaries, 2 = escape at surface only fgets(buf, 32,fid); sscanf(buf, "%lf", &xs); // initial launch point fgets(buf, 32,fid); sscanf(buf, "%lf", &ys); // initial launch point fgets(buf, 32,fid); sscanf(buf, "%lf", &zs); // initial launch point fgets(buf, 32,fid); sscanf(buf, "%lf", &xfocus); // xfocus fgets(buf, 32,fid); sscanf(buf, "%lf", &yfocus); // yfocus fgets(buf, 32,fid); sscanf(buf, "%lf", &zfocus); // zfocus fgets(buf, 32,fid); sscanf(buf, "%lf", &ux0); // ux trajectory fgets(buf, 32,fid); sscanf(buf, "%lf", &uy0); // uy trajectory fgets(buf, 32,fid); sscanf(buf, "%lf", &uz0); // uz trajectory fgets(buf, 32,fid); sscanf(buf, "%lf", &radius); // radius fgets(buf, 32,fid); sscanf(buf, "%lf", &waist); // waist // tissue optical properties fgets(buf, 32,fid); sscanf(buf, "%d", &numberOfTissueTypes); // # of tissue types in tissue list for (tissueIterator=1; tissueIterator<=numberOfTissueTypes; tissueIterator++) { fgets(buf, 32, fid); sscanf(buf, "%lf", &muav[tissueIterator]); // absorption coeff [cm^-1] fgets(buf, 32, fid); sscanf(buf, "%lf", &musv[tissueIterator]); // scattering coeff [cm^-1] fgets(buf, 32, fid); sscanf(buf, "%lf", &gv[tissueIterator]); // anisotropy of scatter [dimensionless] } fclose(fid); if (verbose) { printf("simulationTimeFromFile = %0.2f min\n", simulationTimeFromFile); printf("Nx = %d, xSpacing = %0.4f [cm]\n", Nx, xSpacing); printf("Ny = %d, ySpacing = %0.4f [cm]\n", Ny, ySpacing); printf("Nz = %d, zSpacing = %0.4f [cm]\n", Nz, zSpacing); printf("xs = %0.4f [cm]\n", xs); printf("ys = %0.4f [cm]\n", ys); printf("zs = %0.4f [cm]\n", zs); printf("mcflag = %d:\n", mcflag); if (mcflag == 0) printf("launching uniform flat-field beam\n"); else if (mcflag == 1) printf("launching Gaussian beam\n"); else if (mcflag == 2) printf("launching isotropic point source\n"); else if (mcflag == 4) printf("launching DKFZ custom setup\n"); else printf("launching custom setup..\n"); printf("xfocus = %0.4f [cm]\n", xfocus); printf("yfocus = %0.4f [cm]\n", yfocus); printf("zfocus = %0.2e [cm]\n", zfocus); if (launchflag == 1) { printf("Launchflag ON, so launch the following:\n"); printf("ux0 = %0.4f [cm]\n",ux0); printf("uy0 = %0.4f [cm]\n",uy0); printf("uz0 = %0.4f [cm]\n",uz0); } else { printf("Launchflag OFF, so program calculates launch angles.\n"); printf("radius = %0.4f [cm]\n", radius); printf("waist = %0.4f [cm]\n", waist); } if (boundaryflag == 0) printf("boundaryflag = 0, so no boundaries.\n"); else if (boundaryflag == 1) printf("boundaryflag = 1, so escape at all boundaries.\n"); else if (boundaryflag == 2) printf("boundaryflag = 2, so escape at surface only.\n"); else { printf("improper boundaryflag. quit.\n"); //return nullptr; } printf("# of tissues available, Nt = %d\n",numberOfTissueTypes); for (tissueIterator=1; tissueIterator<=numberOfTissueTypes; tissueIterator++) { printf("muav[%d] = %0.4f [cm^-1]\n",tissueIterator,muav[tissueIterator]); printf("musv[%d] = %0.4f [cm^-1]\n",tissueIterator,musv[tissueIterator]); printf(" gv[%d] = %0.4f [--]\n\n",tissueIterator,gv[tissueIterator]); } } /* IMPORT BINARY TISSUE FILE */ totalNumberOfVoxels = Nx*Ny*Nz; if (verbose) printf("%ld (sizeof totalNumberOfVoxels)\n", totalNumberOfVoxels); tissueType = ( char *)malloc(totalNumberOfVoxels*sizeof(char)); /* tissue structure */ muaVector = ( double *)malloc(totalNumberOfVoxels*sizeof(double)); /* tissue structure */ musVector = ( double *)malloc(totalNumberOfVoxels*sizeof(double)); /* tissue structure */ gVector = ( double *)malloc(totalNumberOfVoxels*sizeof(double)); /* tissue structure */ // read binary file tissueFilename = inputFilename +"_T.bin"; fid = fopen(tissueFilename.c_str(), "rb"); if(fid == nullptr) { std::cout << "Could not open file " << tissueFilename << ". Please check the arguments." << std::endl; exit(EXIT_FAILURE); } fread(tissueType, sizeof(char), totalNumberOfVoxels, fid); fclose(fid); // Show tissue on screen, along central z-axis, by listing tissue type #'s. iy = Ny/2; ix = Nx/2; if (verbose) printf("\n central axial profile of tissue types:\n"); for (iz=0; iz* DetectorVoxelList = new std::vector(); ReturnValues() { 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) { 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, int numThreads); int detector_row = 0; bool interpretAsTime = true; int requestedNumberOfPhotons = 100000; float requestedSimulationTime = 0; // in minutes int concurentThreadsSupported = -1; float yOffset = 0; // in mm bool saveLegacy = false; +mitk::PhotoacousticProbe::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 and header (*_T.bin and *_H.mci)", 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", "d", mitkCommandLineParser::Int, "Detector row (0 = all)", "Detectors = {all rows seperately{-1} | all simultaneously{0} | single row{1, 2, ...}}" ); 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 and starts as many jobs as supported concurrently)."); + 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.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" ] ); 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" ) ) { detector_row = us::any_cast( parsedArgs[ "detector" ] ); } 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::PhotoacousticProbe::New(inputXmlProbeDesign, verbose); + if(!m_PhotoacousticProbe->IsValid()) + { + std::cerr << "Xml File was not valid. Simulation failed." << std::endl; + return EXIT_FAILURE; + } + } if(concurentThreadsSupported == 0 || concurentThreadsSupported == -1) { concurentThreadsSupported = std::thread::hardware_concurrency(); if(concurentThreadsSupported == 0) { printf("Could not determine number of available threads. Launching only one."); concurentThreadsSupported = 1; } } InputValues allInput = InputValues(); allInput.LoadValues(inputFilename, yOffset); std::vector allValues(concurentThreadsSupported); std::thread* threads = new std::thread[concurentThreadsSupported]; for(long i=0; i(simulationTimeElapsed).count() << "sec " << std::endl; /**** SAVE Convert data to relative fluence rate [cm^-2] and save. *****/ int totalNumberOfDetectors = allValues[0].DetectorVoxelList->size(); if (verbose) std::cout << "Total number of detectors = " << totalNumberOfDetectors << std::endl; if (verbose) std::cout << "Allocating memory for result ... "; double* finalTotalFluence; std::vector* detectorFluences = new std::vector(); for(int detectorIndex = 0; detectorIndex < totalNumberOfDetectors; detectorIndex++) { detectorFluences->push_back((double*) malloc(allInput.totalNumberOfVoxels*sizeof(double))); } finalTotalFluence = (double *)malloc(allInput.totalNumberOfVoxels*sizeof(double)); if (verbose) std::cout << "[OK]" << std::endl; if (verbose) std::cout << "Cleaning memory for result ..."; for (int i=0; iat(detectorIndex)[i] = 0; } } if (verbose) std::cout << "[OK]" << std::endl; //Add fluence if (verbose) std::cout << "Calculating resulting fluence ... "; double tdx=0, tdy=0, tdz=0; long long tNphotons=0; for(int t=0; tat(detectorIndex)[voxelNumber] += allValues[t].DetectorVoxelList->at(detectorIndex)->fluenceContribution[voxelNumber]; } } } if (verbose) std::cout << "[OK]" << std::endl; if (verbose) { std::cout << "dx = " << tdx << std::endl; std::cout << "dy = " << tdy << std::endl; std::cout << "dz = " << tdz << std::endl; std::cout << "total number of voxels = " << allInput.totalNumberOfVoxels << std::endl; std::cout << "number of photons = " << (int) tNphotons << std::endl; std::cout << "myname = " << allValues[0].myname << std::endl; } // Normalize deposition (A) to yield fluence rate (F). double temp = tdx*tdy*tdz*tNphotons; for (int i=0; iat(detectorIndex)[i] /= temp*allInput.muaVector[i]; } } if (verbose) std::cout << "Saving Total Fluence 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("y-offset", mitk::PropertyPersistenceInfo::New("y-offset")); mitk::IOUtil::SaveImage(resultImage, outputFilename); if (verbose) std::cout << "[OK]" << std::endl; // Save the binary file std::string outputFolderName = outputFilename.substr(0, outputFilename.find(".nrrd")); #ifdef __linux__ mkdir(outputFolderName.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); if (verbose) std::cout << "Created folder " << outputFolderName << " for linux machine."<< std::endl; #else mkdir(outputFolderName.c_str()); if (verbose) std::cout << "Created folder " << outputFolderName << " for windows machine."<< std::endl; #endif outputFilename = outputFolderName + "/legacyTotalFluence_F.bin"; if (verbose) std::cout << "Saving Legacy Total Fluence to " << outputFilename << std::endl; FILE* fid = fopen(outputFilename.c_str(), "wb"); /* 3D voxel output */ fwrite(finalTotalFluence, sizeof(double), allInput.totalNumberOfVoxels, fid); fclose(fid); if (verbose) std::cout << "... [OK]" << std::endl; if (verbose) std::cout << "Saving Fluence Contributions ..." ; for(int detectorIndex = 0; detectorIndex < totalNumberOfDetectors; detectorIndex++) { std::stringstream detectorname(""); detectorname << allValues[0].DetectorVoxelList->at(detectorIndex)->location.x << "," << allValues[0].DetectorVoxelList->at(detectorIndex)->location.y << "," << allValues[0].DetectorVoxelList->at(detectorIndex)->location.z << "FluenceContribution.bin"; // Save the binary file outputFilename = outputFilename + "/" + allValues[0].myname +"_p" + detectorname.str().c_str(); if (verbose) std::cout << " to " << outputFilename << " ... " << std::endl; FILE* fid1 = fopen(outputFilename.c_str(), "wb"); /* 3D voxel output */ fwrite(detectorFluences->at(detectorIndex), sizeof(double), allInput.totalNumberOfVoxels, fid1); fclose(fid1); } if (verbose) std::cout << "[OK]" << std::endl; exit(EXIT_SUCCESS); } /* end of main */ /* CORE FUNCTION */ void runMonteCarlo(InputValues* inputValues, ReturnValues* returnValue, int thread, int numThreads) { if (verbose) printf("Thread %d: Locking Mutex ...", thread); 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_row == -1) { for (int iz = 0; iz < inputValues->Nz; iz++){ for (int ix = 0; ix < inputValues->Nx/2; ix++){ returnValue->DetectorVoxelList->push_back(new DetectorVoxel(initLocation(ix,inputValues->Ny/2,iz,0),inputValues->totalNumberOfVoxels)); } } } else if(detector_row > 0) { if(detector_row > inputValues->Nz) { std::cout << "Requested detector number not valid. Needs to be >1 and <" << inputValues->Nz << std::endl; exit(EXIT_FAILURE); } for(int ix = 3; ix < inputValues->Nx/2; ix++) { try { returnValue->DetectorVoxelList->push_back(new DetectorVoxel(initLocation(ix,inputValues->Ny/2,detector_row,0),inputValues->totalNumberOfVoxels)); } catch(...) { std::cout << "Out of memory at ix= " << ix << "You will need to make your volume smaller or RAM bigger." << std::endl; exit(EXIT_FAILURE); } } } /**** ======================== MAJOR CYCLE ============================ *****/ std::chrono::system_clock::time_point startTime, finishTime; 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; jtotalNumberOfVoxels;j++) returnValue->totalFluence[j] = 0; // ensure F[] starts empty. /**** RUN Launch N photons, initializing each one before progation. *****/ if (verbose) std::cout << "Performing Monte Carlo Simulation for " << inputFilename << std::endl; int photonsForTimeEstimation = 5000; if(interpretAsTime) { if (verbose) printf("Thread %u: requesting %0.2f min\n", thread, inputValues->simulationTimeFromFile); returnValue->Nphotons = photonsForTimeEstimation; // will be updated to achieve desired run time, time_min. } else { returnValue->Nphotons = requestedNumberOfPhotons / numThreads; if (verbose) std::cout <<"Thread "<< thread <<": Nphotons = "<< (int)returnValue->Nphotons << " for unknown simulation time."< photonsForTimeEstimation) && (photonIterator%(returnValue->Nphotons/100) == 0)) { if (verbose) std::cout <<"Current progress for Thread #"<< thread << ": " << ((double) photonIterator)/returnValue->Nphotons*100 << "%" << std::endl; std::cout <<"" << std::endl; } // At photonsForTimeEstimation th photon, update Nphotons to achieve desired runtime (time_min) if(interpretAsTime) { if (requestedSimulationTime == 0) requestedSimulationTime = inputValues->simulationTimeFromFile; if (photonIterator == 1) startTime = std::chrono::system_clock::now(); if (photonIterator == photonsForTimeEstimation) { finishTime = std::chrono::system_clock::now(); auto elapsed = finishTime - startTime; double simulationTimePerPhoton = std::chrono::duration_cast(elapsed).count() / 60000.0 / (photonsForTimeEstimation - 1.0); returnValue->Nphotons = (long)(requestedSimulationTime / simulationTimePerPhoton); if (verbose) std::cout << "Started Thread " << thread << " with " << returnValue->Nphotons << " photons for an estimated simulation time of " << (requestedSimulationTime)<<"min"<launchflag == 1) // manually set launch + if(m_PhotoacousticProbe.IsNotNull()) { - x = inputValues->xs; - y = inputValues->ys; - z = inputValues->zs; - ux = inputValues->ux0; - uy = inputValues->uy0; - uz = inputValues->uz0; + 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::PhotoacousticLightSource::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 // use mcflag + else { - if (inputValues->mcflag == 0) // uniform beam + /* trajectory */ + if (inputValues->launchflag == 1) // manually set launch { - // 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); + x = inputValues->xs; + y = inputValues->ys; + z = inputValues->zs; + ux = inputValues->ux0; + uy = inputValues->uy0; + uz = inputValues->uz0; } - else if (inputValues->mcflag == 5) // Multispectral DKFZ prototype + else // use mcflag { - // set launch point and width of beam - while ((rnd = returnValue->RandomGen(1,0,NULL)) <= 0.0); + 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; + //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? inputValues->ys+1.5 : inputValues->ys-1.5); - z = 0.1; - ux = 0; + 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); + 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); + //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); + 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); + // 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; + //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? inputValues->ys+0.83 : inputValues->ys-0.83); - z = 0.1; - ux = 0; + 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); + 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); + //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); + 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 + // 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. for(unsigned int j = 0; j < returnValue->DetectorVoxelList->size(); j++) { returnValue->DetectorVoxelList->at(j)->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(unsigned int j = 0; j < returnValue->DetectorVoxelList->size(); j++) { returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->push_back(initLocation(ix, iy, iz, absorb)); if((returnValue->DetectorVoxelList->at(j)->location.x == ix) && ((returnValue->DetectorVoxelList->at(j)->location.y == iy) || (returnValue->DetectorVoxelList->at(j)->location.y-1 == iy)) && (returnValue->DetectorVoxelList->at(j)->location.z == iz)) { for(unsigned int routeIndex = 0; routeIndex < returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->size(); routeIndex++) { i = (long)(returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).z*inputValues->Ny*inputValues->Nx + returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).x*inputValues->Ny + returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).y); returnValue->DetectorVoxelList->at(j)->fluenceContribution[i] += returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).absorb; } returnValue->DetectorVoxelList->at(j)->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 blag==1, i.e., photon inside simulation cube for(unsigned int j = 0; j < returnValue->DetectorVoxelList->size(); j++) { returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->push_back(initLocation(ix, iy, iz, absorb)); if((returnValue->DetectorVoxelList->at(j)->location.x == ix) && ((returnValue->DetectorVoxelList->at(j)->location.y == iy) || (returnValue->DetectorVoxelList->at(j)->location.y-1 == iy)) && (returnValue->DetectorVoxelList->at(j)->location.z == iz)) { for(unsigned int routeIndex = 0; routeIndex < returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->size(); routeIndex++) { i = (long)(returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).z*inputValues->Ny*inputValues->Nx + returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).x*inputValues->Ny + returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).y); returnValue->DetectorVoxelList->at(j)->fluenceContribution[i] += returnValue->DetectorVoxelList->at(j)->recordedPhotonRoute->at(routeIndex).absorb; } returnValue->DetectorVoxelList->at(j)->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 < returnValue->Nphotons); /* end RUN */ if (verbose) printf("------------------------------------------------------\n"); if (verbose) std::cout << "Thread "<< thread <<" is finished." << std::endl; } diff --git a/Modules/PhotoacousticSimulation/Testing/CMakeLists.txt b/Modules/PhotoacousticSimulation/Testing/CMakeLists.txt index 2d5e45bb88..153cd81e2e 100644 --- a/Modules/PhotoacousticSimulation/Testing/CMakeLists.txt +++ b/Modules/PhotoacousticSimulation/Testing/CMakeLists.txt @@ -1,2 +1 @@ MITK_CREATE_MODULE_TESTS() - diff --git a/Modules/PhotoacousticSimulation/Testing/Resources/allsources.xml b/Modules/PhotoacousticSimulation/Testing/Resources/allsources.xml new file mode 100644 index 0000000000..36a213c5f9 --- /dev/null +++ b/Modules/PhotoacousticSimulation/Testing/Resources/allsources.xml @@ -0,0 +1,70 @@ + + + + 1 + + POINT + 0.3 + 0.2 + 0.1 + + + + -0.1 + 0.1 + UNIFORM + + + -0.2 + 0.2 + UNIFORM + + + + + 1 + + RECTANGLE + 0.4 + 0.5 + 0.6 + 0.7 + 0.8 + 0.9 + + + + -0.3 + 0.3 + GAUSSIAN + + + -0.4 + 0.4 + GAUSSIAN + + + + + 1 + + CIRCLE + 0.4 + 0.5 + 0.6 + 0.5 + + + + -0.3 + 0.3 + GAUSSIAN + + + -0.4 + 0.4 + GAUSSIAN + + + + diff --git a/Modules/PhotoacousticSimulation/Testing/Resources/circlesource.xml b/Modules/PhotoacousticSimulation/Testing/Resources/circlesource.xml new file mode 100644 index 0000000000..14ef65dd58 --- /dev/null +++ b/Modules/PhotoacousticSimulation/Testing/Resources/circlesource.xml @@ -0,0 +1,25 @@ + + + + 1 + + CIRCLE + 0.3 + 0.2 + 0.1 + 2 + + + + -0.1 + 0.1 + GAUSSIAN + + + -0.2 + 0.2 + UNIFORM + + + + diff --git a/Modules/PhotoacousticSimulation/Testing/Resources/pointsource.xml b/Modules/PhotoacousticSimulation/Testing/Resources/pointsource.xml new file mode 100644 index 0000000000..0ab73fdb8d --- /dev/null +++ b/Modules/PhotoacousticSimulation/Testing/Resources/pointsource.xml @@ -0,0 +1,24 @@ + + + + 1 + + POINT + 0.3 + 0.2 + 0.1 + + + + -0.1 + 0.1 + UNIFORM + + + -0.2 + 0.2 + UNIFORM + + + + diff --git a/Modules/PhotoacousticSimulation/Testing/Resources/rectanglesource.xml b/Modules/PhotoacousticSimulation/Testing/Resources/rectanglesource.xml new file mode 100644 index 0000000000..7c9a4f2f9e --- /dev/null +++ b/Modules/PhotoacousticSimulation/Testing/Resources/rectanglesource.xml @@ -0,0 +1,27 @@ + + + + 1 + + RECTANGLE + 0.3 + 0.2 + 0.1 + 1 + 2 + 3 + + + + -0.1 + 0.1 + UNIFORM + + + -0.2 + 0.2 + UNIFORM + + + + diff --git a/Modules/PhotoacousticSimulation/Testing/Resources/twopointsources.xml b/Modules/PhotoacousticSimulation/Testing/Resources/twopointsources.xml new file mode 100644 index 0000000000..5054ee51c8 --- /dev/null +++ b/Modules/PhotoacousticSimulation/Testing/Resources/twopointsources.xml @@ -0,0 +1,45 @@ + + + + 1 + + POINT + 0.3 + 0.2 + 0.1 + + + + -0.1 + 0.1 + UNIFORM + + + -0.2 + 0.2 + UNIFORM + + + + + 1 + + POINT + 0.4 + 0.5 + 0.6 + + + + -0.3 + 0.3 + GAUSSIAN + + + -0.4 + 0.4 + GAUSSIAN + + + + diff --git a/Modules/PhotoacousticSimulation/Testing/files.cmake b/Modules/PhotoacousticSimulation/Testing/files.cmake index 876d0bfa8a..615e5513a6 100644 --- a/Modules/PhotoacousticSimulation/Testing/files.cmake +++ b/Modules/PhotoacousticSimulation/Testing/files.cmake @@ -1,24 +1,33 @@ set(MODULE_TESTS # IMPORTANT: If you plan to deactivate / comment out a test please write a bug number to the commented out line of code. # # Example: #mitkMyTest #this test is commented out because of bug 12345 # # It is important that the bug is open and that the test will be activated again before the bug is closed. This assures that # no test is forgotten after it was commented out. If there is no bug for your current problem, please add a new one and # mark it as critical. ################## ON THE FENCE TESTS ################################################# # none ################## DISABLED TESTS ##################################################### # none ################# RUNNING TESTS ####################################################### mitkPhotoacousticTissueGeneratorTest.cpp mitkPhotoacousticVectorTest.cpp mitkPhotoacousticVolumeTest.cpp mitkPhotoacousticVesselTest.cpp mitkPhotoacousticVesselTreeTest.cpp mitkPhotoacousticVesselMeanderStrategyTest.cpp + mitkMcxyzXmlTest.cpp +) + +set(RESOURCE_FILES +pointsource.xml +circlesource.xml +rectanglesource.xml +twopointsources.xml +allsources.xml ) diff --git a/Modules/PhotoacousticSimulation/Testing/mitkMcxyzXmlTest.cpp b/Modules/PhotoacousticSimulation/Testing/mitkMcxyzXmlTest.cpp new file mode 100644 index 0000000000..1acedf6e57 --- /dev/null +++ b/Modules/PhotoacousticSimulation/Testing/mitkMcxyzXmlTest.cpp @@ -0,0 +1,243 @@ +/*=================================================================== + +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 + +// us +#include +#include +#include +#include +#include + +#include +#include + +class mitkMcxyzXmlTestSuite : public mitk::TestFixture +{ + CPPUNIT_TEST_SUITE(mitkMcxyzXmlTestSuite); + MITK_TEST(TestCreatePointSourceProbe); + MITK_TEST(TestCreateCircleSourceProbe); + MITK_TEST(TestCreateRectangleSourceProbe); + MITK_TEST(TestCreateTwoPointSourcesProbe); + MITK_TEST(TestCreateAllSourcesProbe); + MITK_TEST(TestValuesAreInBoundsUniformRectangle); + MITK_TEST(TestValuesAreInBoundsGaussianRectangle); + CPPUNIT_TEST_SUITE_END(); + +private: + + mitk::PhotoacousticLightSource::Pointer m_LightSource; + mitk::PhotoacousticProbe::Pointer m_Probe; + std::string m_XmlProbePointSource; + std::string m_XmlProbeCircleSource; + std::string m_XmlProbeRectangleSource; + std::string m_XmlProbeTwoPointSources; + std::string m_XmlProbeAllSources; + +public: + + void setUp() + { + m_LightSource = mitk::PhotoacousticLightSource::New(); + LoadXmlFile("pointsource.xml", &m_XmlProbePointSource); + LoadXmlFile("circlesource.xml", &m_XmlProbeCircleSource); + LoadXmlFile("rectanglesource.xml", &m_XmlProbeRectangleSource); + LoadXmlFile("twopointsources.xml", &m_XmlProbeTwoPointSources); + LoadXmlFile("allsources.xml", &m_XmlProbeAllSources); + } + + void LoadXmlFile(std::string filename, std::string* lines) + { + us::ModuleResource pointSourceXml = us::GetModuleContext()->GetModule()->GetResource(filename); + std::string line; + if(pointSourceXml.IsValid() && pointSourceXml.IsFile()) + { + us::ModuleResourceStream stream(pointSourceXml); + stream.std::istream::imbue(std::locale("en_GB.UTF-8")); + while(std::getline(stream, line)) + { + *lines = *lines + line + " "; + } + } + else + { + MITK_ERROR << "Xml file was not valid"; + } + } + + void TestCreatePointSourceProbe() + { + m_Probe = mitk::PhotoacousticProbe::New(m_XmlProbePointSource.c_str(), true); + CPPUNIT_ASSERT(true == m_Probe->IsValid()); + } + + void TestCreateCircleSourceProbe() + { + m_Probe = mitk::PhotoacousticProbe::New(m_XmlProbeCircleSource.c_str(), true); + CPPUNIT_ASSERT(true == m_Probe->IsValid()); + } + + void TestCreateRectangleSourceProbe() + { + m_Probe = mitk::PhotoacousticProbe::New(m_XmlProbeRectangleSource.c_str(), true); + CPPUNIT_ASSERT(true == m_Probe->IsValid()); + } + + void TestCreateTwoPointSourcesProbe() + { + m_Probe = mitk::PhotoacousticProbe::New(m_XmlProbeTwoPointSources.c_str(), true); + CPPUNIT_ASSERT(true == m_Probe->IsValid()); + } + + void TestCreateAllSourcesProbe() + { + m_Probe = mitk::PhotoacousticProbe::New(m_XmlProbeAllSources.c_str(), true); + CPPUNIT_ASSERT(true == m_Probe->IsValid()); + } + + void TestValuesAreInBoundsUniformRectangle() + { + int MAXIMUM = 2; + int MINIMUM = -2; + + int ANGLE_MAXIMUM = 1; + int ANGLE_MINIMUM = -1; + + m_LightSource->SetAngleXMode(mitk::PhotoacousticLightSource::DistributionMode::UNIFORM); + m_LightSource->SetAngleYMode(mitk::PhotoacousticLightSource::DistributionMode::UNIFORM); + + m_LightSource->SetAngleXMaximum(ANGLE_MAXIMUM); + m_LightSource->SetAngleXMinimum(ANGLE_MINIMUM); + + m_LightSource->SetAngleYMaximum(ANGLE_MAXIMUM); + m_LightSource->SetAngleYMinimum(ANGLE_MINIMUM); + + m_LightSource->SetSpawnLocationX(MINIMUM); + m_LightSource->SetSpawnLocationXLength(2*MAXIMUM); + + m_LightSource->SetSpawnLocationY(MINIMUM); + m_LightSource->SetSpawnLocationYLength(2*MAXIMUM); + + m_LightSource->SetSpawnLocationZ(MINIMUM); + m_LightSource->SetSpawnLocationZLength(2*MAXIMUM); + + m_LightSource->SetSpawnLocationRadius(MAXIMUM); + + m_LightSource->SetVerbose(false); + + m_LightSource->SetSpawnType(mitk::PhotoacousticLightSource::SpawnType::RECTANGLE); + + std::mt19937 rng; + rng.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock:: + now().time_since_epoch()).count()); + std::uniform_real_distribution<> realDist(0, 1); + + for(int i=0, length=10000; iGetNextPhoton(realDist(rng), realDist(rng), realDist(rng), realDist(rng), + realDist(rng), realDist(rng), realDist(rng)); + + CPPUNIT_ASSERT(result.xAngle>=ANGLE_MINIMUM); + CPPUNIT_ASSERT(result.xAngle<=ANGLE_MAXIMUM); + CPPUNIT_ASSERT(result.yAngle>=ANGLE_MINIMUM); + CPPUNIT_ASSERT(result.yAngle<=ANGLE_MAXIMUM); + CPPUNIT_ASSERT(result.zAngle>=0); + CPPUNIT_ASSERT(result.zAngle<=ANGLE_MAXIMUM); + CPPUNIT_ASSERT(result.xPosition>=MINIMUM); + CPPUNIT_ASSERT(result.xPosition<=MAXIMUM); + CPPUNIT_ASSERT(result.yPosition>=MINIMUM); + CPPUNIT_ASSERT(result.yPosition<=MAXIMUM); + CPPUNIT_ASSERT(result.zPosition>=MINIMUM); + CPPUNIT_ASSERT(result.zPosition<=MAXIMUM); + } + } + + void TestValuesAreInBoundsGaussianRectangle() + { + int MAXIMUM = 2; + int MINIMUM = -2; + + int ANGLE_MAXIMUM = 1; + int ANGLE_MINIMUM = -1; + + m_LightSource->SetAngleXMode(mitk::PhotoacousticLightSource::DistributionMode::GAUSSIAN); + m_LightSource->SetAngleYMode(mitk::PhotoacousticLightSource::DistributionMode::GAUSSIAN); + + m_LightSource->SetAngleXMaximum(ANGLE_MAXIMUM); + m_LightSource->SetAngleXMinimum(ANGLE_MINIMUM); + + m_LightSource->SetAngleYMaximum(ANGLE_MAXIMUM); + m_LightSource->SetAngleYMinimum(ANGLE_MINIMUM); + + m_LightSource->SetSpawnLocationX(MINIMUM); + m_LightSource->SetSpawnLocationXLength(2*MAXIMUM); + + m_LightSource->SetSpawnLocationY(MINIMUM); + m_LightSource->SetSpawnLocationYLength(2*MAXIMUM); + + m_LightSource->SetSpawnLocationZ(MINIMUM); + m_LightSource->SetSpawnLocationZLength(2*MAXIMUM); + + m_LightSource->SetSpawnLocationRadius(MAXIMUM); + + m_LightSource->SetVerbose(false); + + m_LightSource->SetSpawnType(mitk::PhotoacousticLightSource::SpawnType::RECTANGLE); + + std::mt19937 rng; + rng.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock:: + now().time_since_epoch()).count()); + std::uniform_real_distribution<> realDist(0, 1); + + for(int i=0, length=10000; iGetNextPhoton(realDist(rng), realDist(rng), realDist(rng), realDist(rng), + realDist(rng), realDist(rng), realDist(rng)); + + CPPUNIT_ASSERT(result.xAngle>=ANGLE_MINIMUM); + CPPUNIT_ASSERT(result.xAngle<=ANGLE_MAXIMUM); + CPPUNIT_ASSERT(result.yAngle>=ANGLE_MINIMUM); + CPPUNIT_ASSERT(result.yAngle<=ANGLE_MAXIMUM); + CPPUNIT_ASSERT(result.zAngle>=0); + CPPUNIT_ASSERT(result.zAngle<=ANGLE_MAXIMUM); + CPPUNIT_ASSERT(result.xPosition>=MINIMUM); + CPPUNIT_ASSERT(result.xPosition<=MAXIMUM); + CPPUNIT_ASSERT(result.yPosition>=MINIMUM); + CPPUNIT_ASSERT(result.yPosition<=MAXIMUM); + CPPUNIT_ASSERT(result.zPosition>=MINIMUM); + CPPUNIT_ASSERT(result.zPosition<=MAXIMUM); + } + } + + void tearDown() + { + m_XmlProbePointSource = ""; + m_XmlProbeCircleSource = ""; + m_XmlProbeRectangleSource = ""; + m_XmlProbeTwoPointSources = ""; + m_XmlProbeAllSources = ""; + m_Probe = nullptr; + } + +}; + +MITK_TEST_SUITE_REGISTRATION(mitkMcxyzXml) diff --git a/Modules/PhotoacousticSimulation/files.cmake b/Modules/PhotoacousticSimulation/files.cmake index 1151cfaf0e..5d9bf15a57 100644 --- a/Modules/PhotoacousticSimulation/files.cmake +++ b/Modules/PhotoacousticSimulation/files.cmake @@ -1,18 +1,20 @@ set(CPP_FILES Utils/mitkPhotoacousticPropertyCalculator.cpp Utils/mitkPhotoacousticSmartVector.cpp Utils/mitkPhotoacousticStatefulObject.cpp Algorithms/mitkPhotoacousticVolume.cpp Algorithms/mitkPhotoacousticTissueGenerator.cpp Algorithms/mitkPhotoacousticVesselTree.cpp Algorithms/mitkPhotoacousticVessel.cpp Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp Algorithms/mitkPhotoacousticNoiseGenerator.cpp Algorithms/mitkPhotoacoustic3dVolume.cpp Algorithms/mitkPhotoacousticFatTissueGenerator.cpp + ProbeDesign/mitkPhotoacousticProbe.cpp + ProbeDesign/mitkPhotoacousticLightSource.cpp ) set(RESOURCE_FILES photoacoustics128x128.png spectralLIB.dat ) diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacoustic3dVolume.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacoustic3dVolume.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacoustic3dVolume.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacoustic3dVolume.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacoustic3dVolume.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacoustic3dVolume.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacoustic3dVolume.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacoustic3dVolume.h diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticFatTissueGenerator.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticFatTissueGenerator.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticFatTissueGenerator.h diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticNoiseGenerator.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticNoiseGenerator.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticNoiseGenerator.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticNoiseGenerator.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticNoiseGenerator.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticNoiseGenerator.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticNoiseGenerator.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticNoiseGenerator.h diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticTissueGenerator.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticTissueGenerator.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticTissueGenerator.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticTissueGenerator.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticTissueGenerator.h diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVessel.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVessel.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVessel.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVessel.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVessel.h diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselMeanderStrategy.h diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselTree.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselTree.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselTree.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVesselTree.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVesselTree.h diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVolume.cpp b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVolume.cpp rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.cpp diff --git a/Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVolume.h b/Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.h similarity index 100% rename from Modules/PhotoacousticSimulation/Algorithms/mitkPhotoacousticVolume.h rename to Modules/PhotoacousticSimulation/src/Algorithms/mitkPhotoacousticVolume.h diff --git a/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticLightSource.cpp b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticLightSource.cpp new file mode 100644 index 0000000000..13a32e32fb --- /dev/null +++ b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticLightSource.cpp @@ -0,0 +1,416 @@ +/*=================================================================== + +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 "mitkPhotoacousticLightSource.h" +#include "math.h" + +mitk::PhotoacousticLightSource::PhotoacousticLightSource() : + m_IsValid(false) +{ +} + +mitk::PhotoacousticLightSource::PhotoacousticLightSource(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::PhotoacousticLightSource::~PhotoacousticLightSource() +{ + +} + +mitk::PhotoacousticLightSource::TransformResult mitk::PhotoacousticLightSource::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::PhotoacousticLightSource::ParsePhotonDirection(TiXmlElement* element) +{ + TiXmlElement* direction = element->FirstChildElement("PhotonDirection"); + if(direction) + { + ParseAngle(direction, "xAngle"); + } + else + { + if(m_Verbose) + std::cerr << "No \"xAngle\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl; + + m_AngleXMinimum = 0; + m_AngleXMaximum = 0; + m_AngleXMode = DistributionMode::UNIFORM; + } + + direction = element->FirstChildElement("PhotonDirection"); + if(direction) + { + ParseAngle(direction, "yAngle"); + } + else + { + if(m_Verbose) + std::cerr << "No \"yAngle\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl; + + m_AngleYMinimum = 0; + m_AngleYMaximum = 0; + m_AngleYMode = DistributionMode::UNIFORM; + } +} + +void mitk::PhotoacousticLightSource::ParseAngle(TiXmlElement* direction, std::string angle) +{ + double minimum; + double maximum; + DistributionMode mode; + + if(m_Verbose) + std::cout << "Parsing " << angle << std::endl; + TiXmlElement* angleElement = direction->FirstChildElement(angle); + if(angleElement) + { + TiXmlElement* angleMin = angleElement->FirstChildElement("min"); + 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 \"min\" tag in xml. Setting min=0" << std::endl; + minimum = 0; + } + + TiXmlElement* angleMax = angleElement->FirstChildElement("max"); + 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 \"max\" tag in xml. Setting max=0" << std::endl; + maximum = 0; + } + + TiXmlElement* angleMode = angleElement->FirstChildElement("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 \"mode\" tag in xml. Setting mode=UNIFORM" << std::endl; + mode = DistributionMode::UNIFORM; + } + } + else + { + if(m_Verbose) + std::cerr << "No \"angleX\" field in xml. Setting to default (0, 0, UNIFORM)." << std::endl; + + maximum = 0; + minimum = 0; + mode = DistributionMode::UNIFORM; + } + + if(strcmp("xAngle", angle.c_str()) == 0) + { + m_AngleXMinimum = minimum; + m_AngleXMaximum = maximum; + m_AngleXMode = mode; + } + else if(strcmp("yAngle", angle.c_str()) == 0) + { + m_AngleYMinimum = minimum; + m_AngleYMaximum = maximum; + m_AngleYMode = mode; + } +} + +void mitk::PhotoacousticLightSource::ParseEnergy(TiXmlElement* element) +{ + TiXmlElement* energy = element->FirstChildElement("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 \"Energy\" field in xml. Setting Energy=1" << std::endl; + m_Energy = 1.0; + } +} + +void mitk::PhotoacousticLightSource::ParsePhotonSpawnArea(TiXmlElement* element) +{ + TiXmlElement* spawnArea = element->FirstChildElement("PhotonSpawnArea"); + if(spawnArea) + { + TiXmlElement* spawnType = spawnArea->FirstChildElement("SpawnType"); + if(spawnType) + { + std::string spawnTypeText = spawnType->GetText(); + if(strcmp("POINT", spawnTypeText.c_str())==0) + { + m_SpawnType = SpawnType::POINT; + if(m_Verbose) + std::cout << "Setting SpawnType=POINT" << std::endl; + } + else if (strcmp("RECTANGLE", spawnTypeText.c_str())==0) + { + m_SpawnType = SpawnType::RECTANGLE; + if(m_Verbose) + std::cout << "Setting SpawnType=RECTANGLE" << std::endl; + } + else if (strcmp("CIRCLE", spawnTypeText.c_str())==0) + { + m_SpawnType = SpawnType::CIRCLE; + if(m_Verbose) + std::cout << "Setting SpawnType=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 \"SpawnType\" element was not provided for this light source. Light source is not valid." << std::endl; + m_IsValid = false; + } + + TiXmlElement* xLocation = spawnArea->FirstChildElement("x"); + if(xLocation) + { + std::string xLocationText = xLocation->GetText(); + m_SpawnLocationX = std::stod(xLocationText); + if(m_Verbose) + std::cout << "Setting x=" << m_SpawnLocationX; + } + else + { + if(m_Verbose) + std::cerr << "No \"x\" field in xml. Setting x=0" << std::endl; + m_SpawnLocationX = 0; + } + + TiXmlElement* yLocation = spawnArea->FirstChildElement("y"); + if(yLocation) + { + std::string yLocationText = yLocation->GetText(); + m_SpawnLocationY = std::stod(yLocationText); + if(m_Verbose) + std::cout << "Setting y=" << m_SpawnLocationY; + } + else + { + if(m_Verbose) + std::cerr << "No \"y\" field in xml. Setting y=0" << std::endl; + m_SpawnLocationY = 0; + } + + TiXmlElement* zLocation = spawnArea->FirstChildElement("z"); + if(zLocation) + { + std::string zLocationText = zLocation->GetText(); + m_SpawnLocationZ = std::stod(zLocationText); + if(m_Verbose) + std::cout << "Setting z=" << m_SpawnLocationZ; + } + else + { + if(m_Verbose) + std::cerr << "No \"z\" field in xml. Setting z=0.1" << std::endl; + m_SpawnLocationZ = 0.1; + } + + TiXmlElement* rLocation = spawnArea->FirstChildElement("r"); + if(rLocation) + { + std::string rLocationText = rLocation->GetText(); + m_SpawnLocationRadius = std::stod(rLocationText); + if(m_Verbose) + std::cout << "Setting r=" << m_SpawnLocationRadius; + } + else + { + if(m_Verbose) + std::cerr << "No \"radius\" field in xml. Setting radius=0" << std::endl; + m_SpawnLocationRadius = 0; + } + + TiXmlElement* xLength = spawnArea->FirstChildElement("xLength"); + if(xLength) + { + std::string xLengthText = xLength->GetText(); + m_SpawnLocationXLength = std::stod(xLengthText); + if(m_Verbose) + std::cout << "Setting xLength=" << m_SpawnLocationXLength << std::endl; + } + else + { + if(m_Verbose) + std::cerr << "No \"xLength\" field in xml. Setting xLength=0" << std::endl; + m_SpawnLocationXLength = 0; + } + + TiXmlElement* yLength = spawnArea->FirstChildElement("yLength"); + if(yLength) + { + std::string yLengthText = yLength->GetText(); + m_SpawnLocationYLength = std::stod(yLengthText); + if(m_Verbose) + std::cout << "Setting yLength=" << m_SpawnLocationYLength << std::endl; + } + else + { + if(m_Verbose) + std::cerr << "No \"yLength\" field in xml. Setting ylength=0" << std::endl; + m_SpawnLocationYLength = 0; + } + + TiXmlElement* zLength = spawnArea->FirstChildElement("zLength"); + if(zLength) + { + std::string zLengthText = zLength->GetText(); + m_SpawnLocationZLength = std::stod(zLengthText); + if(m_Verbose) + std::cout << "Setting zlength=" << m_SpawnLocationZLength << std::endl; + } + else + { + if(m_Verbose) + std::cerr << "No \"zLength\" field in xml. Setting zlength=0" << std::endl; + m_SpawnLocationZLength = 0; + } + + } + else + m_IsValid = false; +} + +mitk::PhotoacousticLightSource::PhotonInformation mitk::PhotoacousticLightSource::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::PhotoacousticLightSource::IsValid() +{ + return m_IsValid; +} diff --git a/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticLightSource.h b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticLightSource.h new file mode 100644 index 0000000000..2ba2c2e6f2 --- /dev/null +++ b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticLightSource.h @@ -0,0 +1,166 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITKPHOTOACOUSTICLIGHTSOURCE_H +#define MITKPHOTOACOUSTICLIGHTSOURCE_H + +#include + +#include "MitkPhotoacousticSimulationExports.h" + +//Includes for smart pointer usage +#include "mitkCommon.h" +#include "itkObject.h" +#include "itkMacro.h" + +#include + +namespace mitk { + +/** + * @brief The PhotoacousticProbe class + * The representation of a PhotoacousticProbe + */ +class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticLightSource : public itk::Object +{ +public: + + mitkClassMacroItkParent(mitk::PhotoacousticLightSource, itk::Object) + itkFactorylessNewMacro(Self) + mitkNewMacro2Param(Self, TiXmlElement*, bool) + + + enum SpawnType + { + POINT, RECTANGLE, CIRCLE + }; + + enum DistributionMode + { + UNIFORM, GAUSSIAN + }; + + struct PhotonInformation + { + double xPosition; + double yPosition; + double zPosition; + double xAngle; + double yAngle; + double zAngle; + }; + + PhotonInformation GetNextPhoton(double rnd1, double rnd2, double rnd3, + double rnd4, double rnd5, double gau1, double gau2); + + bool IsValid(); + + PhotoacousticLightSource(TiXmlElement* element, bool verbose); + PhotoacousticLightSource(); + virtual ~PhotoacousticLightSource(); + + void ParseAngle(TiXmlElement* direction, std::string angle); + + itkGetMacro(SpawnType, SpawnType) + itkSetMacro(SpawnType, SpawnType) + + itkGetMacro(SpawnLocationX, double) + itkSetMacro(SpawnLocationX, double) + + itkGetMacro(SpawnLocationY, double) + itkSetMacro(SpawnLocationY, double) + + itkGetMacro(SpawnLocationZ, double) + itkSetMacro(SpawnLocationZ, double) + + itkGetMacro(SpawnLocationXLength, double) + itkSetMacro(SpawnLocationXLength, double) + + itkGetMacro(SpawnLocationYLength, double) + itkSetMacro(SpawnLocationYLength, double) + + itkGetMacro(SpawnLocationZLength, double) + itkSetMacro(SpawnLocationZLength, double) + + itkGetMacro(SpawnLocationRadius, double) + itkSetMacro(SpawnLocationRadius, double) + + itkGetMacro(Energy, double) + itkSetMacro(Energy, double) + + itkGetMacro(AngleXMinimum, double) + itkSetMacro(AngleXMinimum, double) + + itkGetMacro(AngleXMaximum, double) + itkSetMacro(AngleXMaximum, double) + + itkGetMacro(AngleYMinimum, double) + itkSetMacro(AngleYMinimum, double) + + itkGetMacro(AngleYMaximum, double) + itkSetMacro(AngleYMaximum, double) + + itkGetMacro(AngleXMode, DistributionMode) + itkSetMacro(AngleXMode, DistributionMode) + + itkGetMacro(AngleYMode, DistributionMode) + itkSetMacro(AngleYMode, DistributionMode) + + itkGetMacro(Verbose, bool) + itkSetMacro(Verbose, bool) + +protected: + + const double TWO_PI = 2.0*3.14159265358979323846; + + SpawnType m_SpawnType; + double m_SpawnLocationX; + double m_SpawnLocationY; + double m_SpawnLocationZ; + double m_SpawnLocationXLength; + double m_SpawnLocationYLength; + double m_SpawnLocationZLength; + double m_SpawnLocationRadius; + + double m_Energy; + + double m_AngleXMinimum; + double m_AngleXMaximum; + double m_AngleYMinimum; + double m_AngleYMaximum; + DistributionMode m_AngleXMode; + DistributionMode m_AngleYMode; + + bool m_IsValid; + bool m_Verbose; + + struct TransformResult + { + double z0; + double z1; + }; + + void ParsePhotonDirection(TiXmlElement* element); + void ParseEnergy(TiXmlElement* element); + void ParsePhotonSpawnArea(TiXmlElement* element); + + TransformResult BoxMuellerTransform(double u1, double u2, double mu, double sigma); + +}; + +} + +#endif // MITKPHOTOACOUSTICLIGHTSOURCE_H diff --git a/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticProbe.cpp b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticProbe.cpp new file mode 100644 index 0000000000..e65b026001 --- /dev/null +++ b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticProbe.cpp @@ -0,0 +1,134 @@ +/*=================================================================== + +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 "mitkPhotoacousticProbe.h" + +mitk::PhotoacousticProbe::PhotoacousticProbe(std::string file, bool verbose) : + m_TotalEnergy(0), + m_Verbose(verbose) +{ + TiXmlDocument xmlDoc(file); + bool success = xmlDoc.LoadFile(TIXML_ENCODING_UTF8); + if(m_Verbose) + { + std::cout << "reading of document was " << (success ? "" : "not ") << "successful" << std::endl; + std::cout << "Content of the Xml file:" << std::endl; + xmlDoc.Print(); + } + if(success) + { + InitProbe(xmlDoc); + } + else + { + if(m_Verbose) + { + std::cerr << "Could not load xml file" << std::endl; + } + } +} + +mitk::PhotoacousticProbe::PhotoacousticProbe(const char* fileStream, bool verbose) : + m_TotalEnergy(0), + m_Verbose(verbose) +{ + TiXmlDocument xmlDoc; + const char* success = xmlDoc.Parse(fileStream, 0, TIXML_ENCODING_UTF8); + if(m_Verbose) + { + std::cout << "reading document was " << (success==nullptr ? "" : "not ")<< "successful (" << (success==nullptr ? "NULL" : success) << ")" << std::endl ; + std::cout << "Content of the Xml file:" << std::endl; + xmlDoc.Print(); + } + if(success == nullptr || atoi(success)==0) + { + InitProbe(xmlDoc); + } + else + { + if(m_Verbose) + { + std::cerr << "Could not load xml file" << std::endl; + } + } +} + +mitk::PhotoacousticProbe::~PhotoacousticProbe() +{ + +} + +mitk::PhotoacousticLightSource::PhotonInformation mitk::PhotoacousticProbe::GetNextPhoton(double rng1, double rnd2, double rnd3, double rnd4, + double rnd5, double rnd6, double rnd7, double rnd8) +{ + rng1 = rng1*m_TotalEnergy; + double currentEnergy = 0; + for(mitk::PhotoacousticLightSource::Pointer lightSource : m_LightSources) + { + currentEnergy += lightSource->GetEnergy(); + if(currentEnergy>=rng1) + return lightSource->GetNextPhoton(rnd2, rnd3, rnd4, rnd5, rnd6, rnd7, rnd8); + } + + //Last resort: If something goes wrong, return a position from the first source. + return m_LightSources[0]->GetNextPhoton(rnd2, rnd3, rnd4, rnd5, rnd6, rnd7, rnd8); +} + +bool mitk::PhotoacousticProbe::IsValid() +{ + return m_IsValid; +} + +void mitk::PhotoacousticProbe::InitProbe(TiXmlDocument xmlDoc) +{ + + m_IsValid = true; + + TiXmlElement* root= xmlDoc.FirstChildElement("Probe"); + if(root) + { + for(TiXmlElement* element = root->FirstChildElement("LightSource"); + element !=nullptr; element = element->NextSiblingElement("LightSource")) + { + mitk::PhotoacousticLightSource::Pointer lightSource = mitk::PhotoacousticLightSource::New(element, m_Verbose); + if(lightSource.IsNotNull() && lightSource->IsValid()) + { + m_LightSources.push_back(lightSource); + m_TotalEnergy += lightSource->GetEnergy(); + } + else + { + m_IsValid = false; + } + } + } + else + { + m_IsValid = false; + } + + if(!m_IsValid) + { + std::cerr << "Creation of a valid Photoacoustic Probe failed." << std::endl; + } + else + { + if(m_Verbose) + { + std::cout << "Successfully created Photoacoustic Probe." << std::endl; + } + } +} diff --git a/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticProbe.h b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticProbe.h new file mode 100644 index 0000000000..03e8571352 --- /dev/null +++ b/Modules/PhotoacousticSimulation/src/ProbeDesign/mitkPhotoacousticProbe.h @@ -0,0 +1,68 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#ifndef MITKPHOTOACOUSTICPROBE_H +#define MITKPHOTOACOUSTICPROBE_H + +#include + +#include + +//Includes for smart pointer usage +#include "mitkCommon.h" +#include "itkLightObject.h" + +#include "mitkPhotoacousticLightSource.h" +#include +#include + +namespace mitk { + +/** + * @brief The PhotoacousticProbe class + * The representation of a PhotoacousticProbe + */ +class MITKPHOTOACOUSTICSIMULATION_EXPORT PhotoacousticProbe : public itk::LightObject +{ +public: + + mitkClassMacroItkParent(mitk::PhotoacousticProbe, itk::LightObject) + mitkNewMacro2Param(Self, std::string, bool) + mitkNewMacro2Param(Self, const char*, bool) + + mitk::PhotoacousticLightSource::PhotonInformation GetNextPhoton(double rng1, double rnd2, double rnd3, double rnd4, + double rnd5, double rnd6, double rnd7, double rnd8); + + bool IsValid(); + + PhotoacousticProbe(std::string xmlFile, bool verbose); + PhotoacousticProbe(const char* fileStream, bool verbose); + virtual ~PhotoacousticProbe(); + +protected: + + std::vector m_LightSources; + bool m_IsValid; + double m_TotalEnergy; + bool m_Verbose; + + void InitProbe(TiXmlDocument document); + +}; + +} + +#endif // MITKPHOTOACOUSTICPROBE_H diff --git a/Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticPropertyCalculator.cpp b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticPropertyCalculator.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticPropertyCalculator.cpp rename to Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticPropertyCalculator.cpp diff --git a/Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticPropertyCalculator.h b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticPropertyCalculator.h similarity index 100% rename from Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticPropertyCalculator.h rename to Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticPropertyCalculator.h diff --git a/Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticSmartVector.cpp b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticSmartVector.cpp rename to Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.cpp diff --git a/Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticSmartVector.h b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.h similarity index 100% rename from Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticSmartVector.h rename to Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticSmartVector.h diff --git a/Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticStatefulObject.cpp b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticStatefulObject.cpp similarity index 100% rename from Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticStatefulObject.cpp rename to Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticStatefulObject.cpp diff --git a/Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticStatefulObject.h b/Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticStatefulObject.h similarity index 100% rename from Modules/PhotoacousticSimulation/Utils/mitkPhotoacousticStatefulObject.h rename to Modules/PhotoacousticSimulation/src/Utils/mitkPhotoacousticStatefulObject.h