montecarlomaplet.maplet

A Monte-Carlo Maplet for the Study of the Optical Properties of Biological Tissues

[Button("Proceed to numerical results",'background'=white,'foreground'=black, 'font'='F1', RunWindow('W1')),
Button("Proceed to graphical results",'background'=white,'foreground'=black, 'font'='F1', RunWindow('W2'))]),
Button['B3']( "Close Window", Shutdown() )]
)),   # end of boxcolumn,boxrow,window

Window['Win']('title'= "Explanation of Input",
BoxRow(
      BoxColumn( 
                 TextBox['TF_I']('editable' = 'false',16..60,  
 " The INPUT WINDOW consists of two columns where data must be entered in order for 
   the Monte Carlo simulation to proceed. 
   On the left column  the entries are:
     number of photons (integer - >500 recommended)
     entry angle (angle in degrees with respect to layer's surface) 
     number of tissue layers (integer)
     layers indices of refraction     (list) 
     layers absorption coefficients   (list)
     layers scattering coefficients   (list)
     layers anisotropy coefficient    (list)
     layers depths                    (list)
  Note that the lists have as many entries as the number of layers; the list format requires that 
  entries be separated by commas and enclosed in square brackets. Ranges of acceptable input values 
  are given in the literature (eg., L.Wang, S.L.Jacques and L.Zheng, Computer Methods and 
                                     Programs in Biomedicine, vol 47, (1995), 131-146)
  
  On the right column, three options, for the layers geometry type, are presented: 
      cylindrical, 
      prismatic and 
      infinite layers.
  If the cylindrical geometry is chosen then the value of the maximum radius has to be given. 
  On the other hand, if the geometry chosen is prismatic, then one needs to give the maximum 
  values of the coordinates x and y. The infinite geometry requires no extra information.
  
  If only the numerical results are required, one should check the appropriate box which 
  speeds up the simulation. Recommended if the number of photons is > 2000.

  Note that the following number of points have been established:
  Na=30 for the exit angle grid with resolution Pi/2Na, Nz=20 with resolution totald/Nz for 
  the depth direction grid (where totald is the sum of the thicknesses of all tissue layers)
  and Nr=300 with resolution dr=0.05 for the radial direction grid.
  
 TO CHANGE THE GRIDS, CLOSE THE MAPLET, GO TO THE PROCEDURE SIMPROP AND CHANGE THE VALUES THERE.

  To start the simulation press the corresponding button.

 CLOSE THIS WINDOW BEFORE PROCEEDING!",   'font'=Font(times,bold,16 ) 
                                 )       # closing TextBox 
))),
Window['W1']('title'= "Numerical Results Window",
BoxRow(
BoxColumn(
         Label("RESULTS OF THE MONTE CARLO SIMULATION",'font'=Font(Arial,bold,16)),

BoxRow('border'='true', 'inset'=0, 
          'font'='F1',                
[MathMLViewer['ML1']('foreground'=black,'background'=white,'height'=25,
          'width'=220),
Button("SPECULAR REFLECTANCE VALUE",'foreground'=black, 'background'=white,'font'='F1',Evaluate('ML1' = 'NumRes1()')) 
 ])
,
BoxRow('border'='true', 'inset'=0, 
           'font'='F1',                
[MathMLViewer['ML2']('foreground'=black,'background'=white,'height'=25,
          'width'=220),
Button("DIFFUSE REFLECTANCE VALUE",'foreground'=black, 'background'=white,'font'='F1',Evaluate('ML2' = 'NumRes2()')) 
])
,
BoxRow('border'='true', 'inset'=0, 
          'font'='F1',                
[MathMLViewer['ML3']('foreground'=black,'background'=white,'height'=25,
          'width'=220),
Button("ABSORBED FRACTION VALUE", 'foreground'=black, 'background'=white,'font'='F1',Evaluate('ML3' = 'NumRes3()')) 
])
,
BoxRow('border'='true', 'inset'=0, 
          'font'='F1',                
[MathMLViewer['ML4']('foreground'=black,'background'=white,'height'=25,
          'width'=220),
Button("TRANSMITTANCE VALUE",'foreground'=black,  'background'=white,'font'='F1',Evaluate('ML4' = 'NumRes4()')), 
Button("CLEAR BOXES", 
Action(
SetOption('ML1' = ""),SetOption('ML2' = ""),SetOption('ML3' = ""),SetOption('ML4' = ""))) 
])

))),


Window['W2']('title'= "Graphical Results Window",
BoxColumn(
  BoxRow('vscroll'='always','hscroll'='always','inset'=0,'spacing'=0,'valign'='top',
     Label("Absorption versus depth Plot      ",'font'=Font(Arial,bold,14)),
 [Plotter['PLA'](Evaluate( 'PLA '= 'GraphMCA()' )),
            Button['B00']("SHOW PLOT",'foreground'=black,  'background'=white, 'font'='F0','enabled'='true',
                                              Action( Evaluate( 'PLA' = 'GraphMCA()'  )) )]
 ),  # end of BoxColumn
BoxRow('vscroll'='always','hscroll'='always','inset'=0,'spacing'=0,'valign'='top',
     Label("Other graphs      ",'font'=Font(Arial,bold,18)),
[ Button("Reflectance Graphs",'foreground'=black,'background'=white, 'font'='F1', RunWindow('W2R')),
Button("Transmittance Graphs",'foreground'=black,'background'=white, 'font'='F1', RunWindow('W2T'))]
))),


Window['W2R']('title'= "Window for Plots for Reflectance ",
BoxRow(
BoxColumn('vscroll'='always','hscroll'='always','inset'=0,'spacing'=0,'valign'='top',
     Label("Reflectance versus radial distance      ",'font'=Font(Arial,bold,16)),
 Plotter['PL2'](Evaluate( 'PL2' = 'GraphMC2()' )),
            Button['B22']("SHOW PLOT",'foreground'=black,  'background'=white, 'font'='F0','enabled'='true',
                                              Action( Evaluate( 'PL2' = 'GraphMC2()'  ) ))
 ),  # end of BoxColumn
BoxColumn('vscroll'='always','hscroll'='always','inset'=0,'spacing'=0,'valign'='top',
      Label("Angular resolution of Reflectance      ",'font'=Font(Arial,bold,16)),
 Plotter['PL4'](Evaluate( 'PL4' = 'GraphMC4()' )),
            Button['B44']("SHOW PLOT",'foreground'=black,  'background'=white, 'font'='F0',  'enabled'='true',
                                              Action(Evaluate( 'PL4' = 'GraphMC4()'  ) ))
 )  # end of BoxColumn
)),
Window['W2T']('title'= "Window for Plots for Transmittance",
BoxRow(
BoxColumn('vscroll'='always','hscroll'='always','inset'=0,'spacing'=0,'valign'='top',
            Label("Transmittance versus radial distance      ",'font'=Font(Arial,bold,16)),
 Plotter['PL1'](Evaluate( 'PL1' = 'GraphMC1()' )),
            Button['B11']("SHOW PLOT",'foreground'=black,  'background'=white, 'font'='F0', 'enabled'='true',
                                              Action( Evaluate( 'PL1' = 'GraphMC1()'  )) )
 ),  # end of BoxColumn
BoxColumn('vscroll'='always','hscroll'='always','inset'=0,'spacing'=0,'valign'='top',
            Label("Angular Resolution of Transmittance      ",'font'=Font(Arial,bold,16)),
 Plotter['PL3'](Evaluate( 'PL3' = 'GraphMC3()' )),
            Button['B33']("SHOW PLOT",'foreground'=black,  'background'=white, 'font'='F0',  'enabled'='true',
                                              Action( Evaluate( 'PL3' = 'GraphMC3()'  ) ))
 )  # end of BoxColumn
)),

Window['Wout']('title'= "Explanation of Output",
BoxRow(
BoxColumn( 
                 TextBox['TF_O']('editable' = 'false',20..65,  
 "     The output of the simulation can be obtained in numerical or graphical format.
    Activating the button [Proceed to numerical results] opens a new window where the values of the  
          Specular Reflectance        (TR) 
          Diffuse Reflectance         (RD)
          Absorbed Fraction and       (AA)
          Transmittance               (TD)
    retrieved from the Maple internal file  MCnumerical.m   are displayed.
        
       On the other hand, the button [Proceed to graphical results] opens a window where the 
     graph of absorption versus layers depth (stored in MCgraphA.m) can be seen.
     At the bttom of this same window, there are two buttons which lead to the display of graphs 
     of the radial-resolved and angular-resolved reflectance and transmittance. 
  
        These radial resolved and angular resolved graphs are retrieved from the Maple internal files
     MCgraphRad.m   and  MCgraphAng.m  respectively.   

     Furthermore, in case the points of these graphs are needed they can be retrieved from the files:
          Abs_graph_points.m       --- where the points for the absorption versus depth graph points are stored
          Ref_ang_graph_points.m   --- where the points for the reflectance versus azimuthal angle graph points are stored
          Tra_ang_graph_points.m   --- where the points for the tranmitance versus azimuthal angle graph points are stored
          Ref_rad_graph_points.m   --- where the points for the reflectance versus radial coordinates graph points are stored
          Tra_rad_graph_points.m   --- where the points for the tranmitance versus radial coordinates graph points are stored   

  CLOSE THIS WINDOW BEFORE PROCEEDING!  ",   'font'=Font(times,bold,14)
                                ) # closing TextBox
)))
):  #  closing Maplet
end use:

# 2. Procedures for simulation when running the Maplet
SimulationMC:=proc()
global layers, photons,geom,n, mua,mus,g,d,totald,rmax,xymax,xmax,ymax,entryang;
local i,nn,mmua,mmus,gg,dd,g1,g2,g3,g4,graphs;

use Maplets[Tools] in 
g1:=Get('ChB1'):
g2:=Get('ChB2'):
g3:=Get('ChB3'):
g4:=Get('ChB4'):
if g4=true then graphs:=0; else graphs:=1;fi;
if g1=true then geom:=cylindrical;
rmax:=Get('TF_r'::integer);
elif g2=true then geom:=prismatic; 
xymax:=Get('TF_xy'::list);
xmax:=xymax[1];
ymax:=xymax[2];
else  geom:=infinite;fi;
entryang:=Get('TF_0'::anything);
layers:=Get('TF_1'::integer);
photons:=Get('TF_2':: integer);
nn:=Get ('TF_3'::list);
mmus:=Get ('TF_4'::list);
mmua:=Get ('TF_5'::list);
gg:=Get ('TF_6'::list);
dd:=Get ('TF_7'::list);
end use:
                            # Initialize layer property arrays
 n := array(0..layers+1):
 mua := array(1..layers):   # changed the first index of the array from 0 to 1
 mus := array(1..layers):   # same change as above
 g := array(1..layers):
 d := array(0..layers): 
totald:=0;
for i to layers do
n[i]:=nn[i];
mua[i]:=mmua[i];
mus[i]:=mmus[i];
g[i]:=gg[i];
d[i]:=dd[i];
totald:=totald+d[i];
od;
MonteCarloInMaple(graphs, geom,photons,entryang):
RETURN("Simulation finished - Proceed to Results "); 
end proc:
NumRes1:=proc()
read `MCresults.m`:
Digits:=4;RETURN(TR/1.0);end:
NumRes2:=proc()
read `MCresults.m`:
Digits:=4;RETURN(RD/1.0);end:
NumRes3:=proc()
read `MCresults.m`:
Digits:=4;RETURN(AA/1.0);end:
NumRes4:=proc()
read `MCresults.m`:
Digits:=4;RETURN(TD/1.0);end:
GraphMCA:=proc()
read `MCgraphA.m`:
RETURN(display(PPa));end:
GraphMC1:=proc()
read `MCgraphRad.m`:
RETURN(display(PPrT));end:
GraphMC2:=proc()
read `MCgraphRad.m`:
RETURN(display(PPrR));end:
GraphMC3:=proc()
read `MCgraphAng.m`:
RETURN(display(PPaT));end:
GraphMC4:=proc()
read `MCgraphAng.m`:
RETURN(display[plots](PPaR));end:
# III - Running the simulation 
# 1. Running the simulation from the Maplet (it rquires removing the symbol # appearing at the beginning of the next command line)
Maplets[Display](MonteCarloMaplet);
# 2. Running the simulation directly from the worksheet. Check/change the numerical values (purple statements) inside the following procedure.
# 
RUN_WORKSHEET_SIMULATION:=proc()
global layers, photons,geom,n,mua,mus,g,d,totald,rmax,xymax,xmax,ymax,entryang;
local i,nn,mmua,mmus,gg,dd,g1,g2,g3,g4,graphs,geo;

photons:=1000;    # number of photons  - integer
entryang:=90;   # angle the photon makes with the surface layer - real number between  0 and 90o
layers:=2;      # Number of the tissue layers - integer
nn:=[1.5,1.33];     # indices of refraction of the layers - list
mmua:=[0.1,0.3];    # coefficient of absorption of the layers - list (entries are asumed in 1/cm)
mmus:=[10,20];     # coefficient of scattering of the layers - list (entries are assumed in 1/cm)
gg:=[0.93,0.87];     # anisotropy of the layers - list
dd:=[0.1,0.2];      # depth of the layers - list (entries are assumed in cm)
      geo:=[infinite,prismatic,cylindrical];  # choice of three types of layers' geometry
geom:=geo[1];   # here made choice 1.
     if geom=geo[3] then rmax:=1;
     elif geom=geo[2] then xmax:=16; ymax:=16;
     fi;        # give dimensions of layers
graphs:=1;      
# enter graphs:=0 for computing only numerical results; graphs:=1 for numerical and graphical 
# results
 n := array(0..layers+1):
 mua := array(1..layers):   # changed the first index of the array from 0 to 1
 mus := array(1..layers):   # same change as above
 g := array(1..layers):
 d := array(0..layers): 
totald:=0;
for i to layers do
n[i]:=nn[i];
mua[i]:=mmua[i];
mus[i]:=mmus[i];
g[i]:=gg[i];
d[i]:=dd[i];
totald:=totald+d[i];
od;
MonteCarloInMaple(graphs,geom,photons,entryang);
end:
# remove the symbol # appearing at the beginning of the next coomand line.
#RUN_WORKSHEET_SIMULATION();
# IV - Retrieving stored data 
# 1. Reading in the files 
read `MCresults.m`:   # where numerical results are stored
read `MCgraphA.m`:    # where the depth resolved absorption graph is stored
read `MCgraphAng.m`:  # where the angular resolved reflectance/transmittance graphs are stored
read `MCgraphRad.m`:  # where the radial resolved reflectance/transmittance graphs are stored
read `Abs_graph_points.m`:# where the absorption versus dept plot points are stored
read`Tra_ang_graph_points.m`:# where the transmittance versus angle plot points are stored
read`Tra_rad_graph_points.m`:# where the transmittance versus radius plot points are stored
read`Ref_ang_graph_points.m`:# where the reflectance versus angle plot points are stored
read`Ref_rad_graph_points.m`:# where the reflectance versus radius plot points are stored


# 2. Examples
# 2.1 Numerical results
Digits:=4:
print("Specular Reflectance =",TR/1.0);
print("Diffuse Reflectance =",RD/1.0);
print("Absorbed fraction =" ,AA/1.0);
print("Transmittance= ",TD/1.0);

# 2.2 Plotting from the stored point coordinates
pointplot(ppaR, labels=["Angle (rad)", Reflectance], axes=boxed);

pointplot(pprR, labels=["radius(cm)", Reflectance], axes=boxed);
pointplot(ppa, labels=["depth(cm)", Absorption], axes=boxed);

pointplot(ppaT, labels=["Angle (rad)", Transmittance], axes=boxed);
pointplot(pprT, labels=["radius(cm)", Transmittance], axes=boxed);
# 2.3 Displaying the stored plots
PPa;
PPaR;
PPrR;
PPaT;
PPrT;