Commit: e98db29ddb5c4dc67104c0e4c5adfb3aff342110
Author: Randy Palamar
Date: Tue, 26 Mar 2024 11:12:00 -0600
init
Diffstat:
| A | Photon.cpp | | | 248 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | Photon.h | | | 44 | ++++++++++++++++++++++++++++++++++++++++++++ |
| A | Utilities.cpp | | | 94 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | Utilities.h | | | 140 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| A | main.cpp | | | 143 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
5 files changed, 669 insertions(+), 0 deletions(-)
diff --git a/Photon.cpp b/Photon.cpp
@@ -0,0 +1,248 @@
+// Photon.cpp: implementation of the Photon class.
+//
+//////////////////////////////////////////////////////////////////////
+
+#include "Photon.h"
+#include "Utilities.h"
+
+#include <math.h>
+/////////////////////////////////////////
+
+
+double** Photon::Rd_xy=NULL;
+double Photon::Rsp=0.0;
+/////////////////////////////////////////////////////////////////////////
+
+
+//////////////////////////////////////////////////////////////////////
+// Construction/Destruction
+//////////////////////////////////////////////////////////////////////
+
+Photon::Photon()
+{
+ w=1.0;
+ nScatter=0;
+ dead=false;
+ posx=0.0;
+ posy=0.0;
+ posz=0.0;
+////////incident direction
+ ux=sin(th_in)*n0/n;
+ uy=0.0;
+ uz=sqrt(1-ux*ux);
+//////////////////////////////////////
+ s=0.0;
+ db=0.0;
+}
+
+Photon::~Photon()
+{
+
+}
+
+void Photon::Launch()
+{
+ double rsp, sinthi, sintht, costhi, costht, raux1, raux2;
+ sinthi=sin(th_in);
+ costhi=cos(th_in);
+ sintht=ux;
+ costht=uz;
+ raux1=sinthi*costht-costhi*sintht;
+ raux1=raux1*raux1;
+ raux2=sinthi*costht+costhi*sintht;
+ raux2=raux2*raux2;
+ rsp=0.5*(raux1/raux2+raux1*(1-raux2)/raux2/(1-raux1));
+ Rsp+=w*rsp;
+ w=w*(1-rsp);
+}
+
+void Photon::CalStep()
+{
+ double ran;
+ if (s<ZERO)
+ {
+ do
+ {
+ ran=RandomNum();
+ } while(ran<=0||ran>1);
+ s=-log(ran); //dimensionless,
+ }
+}
+
+void Photon::Caldb()
+{
+ if (uz<-ZERO) //travel up to reach upper layer
+ {
+ db=-posz/uz;
+ }
+ else if (uz>ZERO) //travel down to reach bottom layer
+ {
+ db=(d-posz)/uz;
+ }
+ else
+ {
+ db=BIGVAL; // A big value
+ }
+}
+
+bool Photon::HitBoundary()
+{
+ if(db*mut<=s)
+ return true;
+ else
+ return false;
+}
+
+void Photon::Move(double ds) //ds has dimension in cm
+{
+ posx+=ds*ux;
+ posy+=ds*uy;
+ posz+=ds*uz;
+ s=s-ds*mut;
+}
+
+void Photon::Absorb() // calculate absorption
+{
+ double dw;
+ dw=w*mua/mut;//3.24
+ w-=dw;
+ if (w<0)
+ {
+ dead=true;
+ w=0.0;
+ dw=w+dw; //dw equals the small weight dw=(w-dw)+dw
+ }
+}
+
+void Photon::Scatter() // calculate next propagation direction
+{
+ if (!dead)
+ {
+ double costheta,fei,ran;
+ do
+ {
+ ran=RandomNum();
+ } while(ran<0||ran>=1);
+ if (g!=0)
+ {
+ double aa=(1-g*g)/(1-g+2*g*ran);
+ costheta=(1+g*g-aa*aa)/(2*g);
+ if(costheta<-1) costheta=-1;
+ else if (costheta>1) costheta=1;
+ }
+ else
+ {
+ costheta=2*ran-1;
+ }//(3.28)
+ double sintheta,sinfei,cosfei;
+ sintheta=sqrt(1-costheta*costheta);
+ do
+ {
+ ran=RandomNum();
+ } while(ran<0||ran>=1);
+ fei=2*PI*ran;//(3.29)
+ cosfei=cos(fei);
+ if (fei<PI)
+ {
+ sinfei=sqrt(1-cosfei*cosfei);
+ }
+ else
+ {
+ sinfei=-sqrt(1-cosfei*cosfei);
+ }
+ double uux,uuy,uuz;
+ if ((SIGN(uz)*uz)<=(1-ZERO))
+ {
+ uux=sintheta*(ux*uz*cosfei-uy*sinfei)/sqrt(1-uz*uz)+ux*costheta;
+ uuy=sintheta*(uy*uz*cosfei+ux*sinfei)/sqrt(1-uz*uz)+uy*costheta;
+ uuz=-sintheta*cosfei*sqrt(1-uz*uz)+uz*costheta;
+ }//(3.24)
+ else //close to normal propagation
+ {
+ uux=sintheta*cosfei;
+ uuy=sintheta*sinfei;
+ uuz=SIGN(uz)*costheta;
+ }//(3.30)
+ ux=uux;
+ uy=uuy;
+ uz=uuz;
+ nScatter+=1;
+ }
+}
+
+void Photon::ReflectTransmit() // deal with photon across boundary
+{
+ double sinalphai,sinalphat,ralphai;
+ sinalphai=sqrt(1-(uz)*(uz));
+ sinalphat=sinalphai*n/n0; //3.35
+ if (sinalphat<ZERO) // close to normal
+ {
+ ralphai=(n-n0)/(n+n0);
+ ralphai=ralphai*ralphai;
+ }
+ else if (sinalphat>1)// Total reflection
+ {
+ ralphai=1;
+ }
+ else
+ {
+ double A,B;
+ A=sinalphai*sqrt(1-sinalphat*sinalphat)-sinalphat*sqrt(1-sinalphai*sinalphai);
+ A=A*A;
+ B=sinalphai*sqrt(1-sinalphat*sinalphat)+sinalphat*sqrt(1-sinalphai*sinalphai);
+ B=B*B;
+ ralphai=0.5*(2*A-A*A-A*B)/(B*(1-A));//(3.36)
+ }
+ double ran;
+ do
+ {
+ ran=RandomNum();
+ }
+ while(ran<0||ran>=1);
+ if (ran<=ralphai) //rebound to current layer
+ {
+ uz=-uz;
+ }
+ else //cross the boundary
+ {
+ if (uz<-ZERO) //hit upper doundary
+ {
+ if(posx>-x_offset&&posx<x_offset&&posy>-y_offset&&posy<y_offset)
+ {
+ Rd_xy[int((posx+x_offset)/dx)][int((posy+y_offset)/dy)]+=w;
+ }
+ dead=true;
+ }
+ else //hit bottom boundary
+ {
+ dead=true;
+ }
+ }
+}
+
+void Photon::Termination() //Russian roulette survial test
+{
+ if (!dead)
+ {
+ double wth=1e-4;
+ double m=10;
+ if (w<=wth)
+ {
+ double e;
+ do
+ {
+ e=RandomNum();
+ }while(e<0||e>1);
+ if (e>(1/m))
+ {
+ w=0;
+ dead=true;
+ }
+ else
+ {
+ w=m*w;
+ }
+ }
+ }
+}
+
diff --git a/Photon.h b/Photon.h
@@ -0,0 +1,44 @@
+// Photon.h: interface for the Photon class.
+// describe the behaviors of photons in Tissue
+//////////////////////////////////////////////////////////////////////
+
+#if !defined(AFX_PHOTON_H__3C1D58D8_26F1_4E3A_9F95_D5051F1FDE29__INCLUDED_)
+#define AFX_PHOTON_H__3C1D58D8_26F1_4E3A_9F95_D5051F1FDE29__INCLUDED_
+
+#if _MSC_VER > 1000
+#pragma once
+#endif // _MSC_VER > 1000
+
+
+
+
+
+class Photon
+{
+public:
+ double posx, posy, posz; //current position of photon
+ double ux, uy, uz; //current direction of photon
+ int nScatter; //# of scattering events experienced
+ double w; //weight
+ double s; //s is the remaining dimensionless step-size
+ double db; //the distance from boundary
+ bool dead; //True if photon is dead, otherwise False
+
+ static double** Rd_xy;
+ static double Rsp; //specular reflectance
+
+
+ Photon();
+ virtual ~Photon();
+ void Launch(); // calculate specular reflectance at the surface
+ void CalStep(); // calculate dimensionless step length
+ void Caldb(); // calculate the travel length to reach the boundary
+ bool HitBoundary(); // True if hitboudary
+ void Move(double ds); //move photon towards boundary
+ void Absorb(); // calculate absorption
+ void Scatter(); // calculate next propagation direction
+ void ReflectTransmit(); // deal with photon across boundary
+ void Termination(); //Russian roulette survial test
+};
+
+#endif // !defined(AFX_PHOTON_H__3C1D58D8_26F1_4E3A_9F95_D5051F1FDE29__INCLUDED_)
diff --git a/Utilities.cpp b/Utilities.cpp
@@ -0,0 +1,94 @@
+
+#include "Utilities.h"
+
+/***********************************************************
+* A random number generator from Numerical Recipes in C. ****/
+#define MBIG 1000000000
+#define MSEED 161803398
+#define MZ 0
+#define FAC 1.0E-9
+float ran3(int &idum)
+{
+ static int inext,inextp;
+ static long ma[56];
+ static int iff=0;
+ long mj,mk;
+ int i,ii,k;
+ if (idum < 0 || iff == 0)
+ {
+ iff=1;
+ mj=MSEED-(idum < 0 ? -idum : idum);
+ mj %= MBIG;
+ ma[55]=mj;
+ mk=1;
+ for (i=1;i<=54;i++)
+ {
+ ii=(21*i) % 55;
+ ma[ii]=mk;
+ mk=mj-mk;
+ if (mk < MZ) mk += MBIG;
+ mj=ma[ii];
+ }
+ for (k=1;k<=4;k++)
+ for (i=1;i<=55;i++)
+ {
+ ma[i] -= ma[1+(i+30) % 55];
+ if (ma[i] < MZ) ma[i] += MBIG;
+ }
+ inext=0;
+ inextp=31;
+ idum=1;
+ }
+ if (++inext == 56)
+ inext=1;
+ if (++inextp == 56)
+ inextp=1;
+ mj=ma[inext]-ma[inextp];
+ if (mj < MZ) mj += MBIG;
+ ma[inext]=mj;
+ return mj*FAC;
+}
+#undef MBIG
+#undef MSEED
+#undef MZ
+#undef FAC
+
+
+/***********************************************************
+ * Generate a random number between 0 and 1. Take a
+ * number as seed the first time entering the function.
+ * The seed is limited to 1<<15.
+ * Wang et al found that when idum is too large, ran3 may return
+ * numbers beyond 0 and 1.
+ ****/
+//#define STANDARDTEST 0 /* testing program using fixed rnd seed. */
+double RandomNum(void)
+{
+ static bool first_time=1;
+ static int idum; /* seed for ran3. */
+ if(first_time)
+ {
+ #if STANDARDTEST /* Use fixed seed to test the program. */
+ idum = - 1;
+ #else
+ idum = -(int)time(NULL)%(1<<15); /* use 16-bit integer as the seed. */
+ #endif
+ ran3(idum);
+ first_time = 0;
+ idum = 1;
+ }
+ return( (double)ran3(idum) );
+}
+
+/////////////////////////////////////////////
+int min(int d1,int d2)
+{
+ if (d1<d2)
+ {
+ return d1;
+ }
+ else
+ {
+ return d2;
+ }
+}
diff --git a/Utilities.h b/Utilities.h
@@ -0,0 +1,140 @@
+/* Providing global variables and functions*/
+#include <math.h>
+#include <new>
+#include <time.h>
+
+//global constants
+#ifndef PI
+ #define PI 3.14159265
+#endif
+
+#ifndef ZERO
+ #define ZERO 1.0e-6 //for small deflection angle
+#endif
+
+#ifndef n0
+ #define n0 1.0 //refractive index of ambient medium
+#endif
+
+#ifndef BIGVAL
+ #define BIGVAL 1e6//big value
+#endif
+
+//global variables
+#ifdef GLOBALORIGIN
+ #define GLOBAL
+#else
+ #define GLOBAL extern
+#endif
+GLOBAL double dx;
+GLOBAL double dy;
+GLOBAL int Nx;
+GLOBAL int Ny;
+GLOBAL double zc;
+GLOBAL double n;
+GLOBAL double mua;
+GLOBAL double mus;
+GLOBAL double mut;
+GLOBAL double g;
+GLOBAL double d;
+
+
+GLOBAL double th_in; // incident angle in degree
+
+GLOBAL double x_offset;
+GLOBAL double y_offset;
+
+//generate random data
+float ran3(int *idum);
+double RandomNum(void);
+
+
+//other auxilliary functions
+int min(int d1, int d2);
+#define SIGN(x) ((x)>=0?1:-1)
+
+
+////////allocte dynamic array/////////////////////////
+using namespace std;
+template <class T>
+T* alloc1D(int dim1, T t)
+{
+ T* p;
+ try
+ {
+ p=new T[dim1];
+ }
+ catch(bad_alloc)
+ {
+ cerr<<"Fail to allocate 1D array!"<<endl;
+ exit(1);
+ }
+ for (int tmp1=0;tmp1<dim1;tmp1++)
+ {
+ p[tmp1]=t;
+ }
+ return p;
+}
+template <class T>
+void dealloc1D(T* p)
+{
+ delete [] p;
+}
+///////////
+template <class T>
+T** alloc2D(int dim1,int dim2,T t)
+{
+ T** p;
+ try
+ {
+ p=new T*[dim1];
+ for (int tmp1=0;tmp1<dim1;tmp1++)
+ {
+ p[tmp1]=alloc1D(dim2,t);
+ }
+ }
+ catch (bad_alloc)
+ {
+ cerr<<"Fail to allocate 2D array!";
+ }
+ return p;
+}
+template <class T>
+void dealloc2D(int dim1,T** p)
+{
+ for (int tmp1=0;tmp1<dim1;tmp1++)
+ {
+ dealloc1D(p[tmp1]);
+ }
+ delete [] p;
+}
+/////////
+template <class T>
+T*** alloc3D(int dim1,int dim2,int dim3,T t)
+{
+ T*** p;
+ try
+ {
+ p=new T**[dim1];
+ for (int tmp1=0;tmp1<dim1;tmp1++)
+ {
+ p[tmp1]=alloc2D(dim2,dim3,t);
+ }
+ }
+ catch (bad_alloc)
+ {
+ cerr<<"Fail to allocate 3D array!";
+ }
+ return p;
+}
+template <class T>
+void dealloc3D(int dim1,int dim2,T*** p)
+{
+ for (int tmp1=0;tmp1<dim1;tmp1++)
+ {
+ dealloc2D(dim2,p[tmp1]);
+ }
+ delete [] p;
+}
+////////////////////////////////////////////////
+
diff --git a/main.cpp b/main.cpp
@@ -0,0 +1,143 @@
+ /* Hybrid-MC by Li Li, Last modified 03/07
+refrence: Wang, "Biomedical optics", chpater3&6 */
+
+
+/* Programming notes
+1. Use Cartesian coordinates.
+2. Incident beam in x-z plane.
+*/
+
+
+#include <iostream>
+#include <fstream>
+#include <string>
+
+#define GLOBALORIGIN
+#include "Utilities.h"
+#undef GLOABLORIGIN
+
+#include "Photon.h"
+
+
+////////////////////////////////////////////////
+void main()
+{
+////////Input Interface////////////////////
+ cout<<"*****************************************"<<endl;
+ cout<<"Hybrid Monte Carlo for Homogeneous Slab with Oblique incidence"<<endl;
+ cout<<"By Li Li, 04/2007"<<endl;
+ cout<<"Refrence: Wang, Biomedical optics, chpater3&6"<<endl;
+ cout<<"*****************************************"<<endl;
+ cout<<endl;
+
+ cout<<"File name(Input:*.mci; Output:*.mco): ";
+ string sFilename, sInfilename, sOutfilename;
+ cin>>sFilename;
+ sInfilename=sFilename+".mci";
+ sOutfilename=sFilename+".mco";
+ ifstream Infile;
+ Infile.open(sInfilename.c_str());
+ if(Infile.fail())
+ {
+ cerr<<"Error! Can't open input file. Err Code 1"<<endl;
+ exit(1);
+ }
+ // double P0, R; //Incident energy/power P0, Beam diameter R
+ // Infile>>P0>>R;
+ int nNum_photons; //# of photons used
+ Infile>>nNum_photons;
+ Infile>>th_in;
+ th_in=th_in/180*PI;
+ Infile>>dx>>dy>>Nx>>Ny;
+ x_offset=dx*Nx/2;
+ y_offset=dy*Ny/2;
+ Infile>>n>>mua>>mus>>g>>d;
+ mut=mua+mus;
+ Infile.close();
+/////////////////////////////////////////////
+
+////////for timing///////////////////////////
+ time_t start,end; //simulation beginning and finishing time
+ time(&start);
+/////////////////////////////////////////////
+ int tmp1, tmp2;
+
+////////initiate recording matrix////////////
+
+ double** Rd_xy=alloc2D(Nx,Ny,0.0);
+
+ Photon::Rd_xy=Rd_xy;
+
+//////////Propagate photons////////////////////////////////
+ for (tmp1=1;tmp1<=nNum_photons;tmp1++)
+ {
+ Photon photon;
+ photon.Launch();
+ do
+ {
+ photon.CalStep();
+ photon.Caldb();
+ if(photon.HitBoundary())
+ {
+ photon.Move(photon.db);
+ photon.ReflectTransmit();
+ }
+ else
+ {
+ photon.Move(photon.s/mut);
+ photon.Absorb();
+ photon.Scatter();
+ }
+ photon.Termination();
+ } while(!photon.dead);
+ if (tmp1%(nNum_photons/10)==0)
+ {
+ cout<<tmp1<<" photons are done!"<<endl;
+ }
+ }
+
+ time(&end);
+ int CompuTime;
+ CompuTime=end-start;
+ cout<<"This simulation take "<<CompuTime<<" seconds!"<<endl;
+/////////////////////////////////////////////
+
+
+////////Write data before convolution////////
+ ofstream f_out(sOutfilename.c_str());
+ f_out<<"/////////////Grid coordinates///////////"<<endl;
+ f_out<<"x (cm): ";
+ for (tmp1=0;tmp1<Nx;tmp1++)
+ {
+ f_out<<(tmp1+0.5)*dx-x_offset<<' ';
+ }
+ f_out<<endl;
+ f_out<<"y (cm): ";
+ for (tmp1=0;tmp1<Ny;tmp1++)
+ {
+ f_out<<(tmp1+0.5)*dy-y_offset<<' ';
+ }
+ f_out<<endl;
+ f_out<<"/////////////////////////////////////////"<<endl;
+
+
+ f_out<<"///////////////Rd_xy/////////////////////"<<endl;
+ double scale=nNum_photons*dx*dy;
+ for (tmp1=0; tmp1<Nx; tmp1++)
+ {
+ for (tmp2=0; tmp2<Ny; tmp2++)
+ {
+ f_out<<Rd_xy[tmp1][tmp2]/scale<<' ';
+ }
+ f_out<<endl;
+ }
+ f_out<<endl;
+ f_out<<"/////////////////////////////////////////"<<endl;
+
+
+
+
+////////deallocate recording matrix//////////
+ dealloc2D(Nx,Rd_xy);
+/////////////////////////////////////////////
+}