GOALS: Gain familiarity with the error characteristics of Monte Carlo simulations

Open VTS_MATLAB folder and click on file short_course_monte_carlo_lab.m.

I. Dependence of Fluence Predictions on Number of Photons Simulated

Goal: This exercise explores how fluence estimates change with the number of photons simulated.

1. Go to the section Example 1a. This example sets up a Monte Carlo simulation for a system with a collimated source impinging on a slab of tissue. The tissue optical properties are specified in the "tissueInput.LayerRegions" struct under "RegionOP".

   tissueInput.LayerRegions = struct(... 
       'ZRange', ... 
       {... 
           [-Inf, 0], ... % air "z" range 
           [0, 100], ... % tissue "z" range 
           [100, +Inf] ... % air "z" range 
       }, ... 
       'RegionOP', ...
       {... 
           [0.0, 1e-10, 1.0, 1.0], ... % air optical properties 
           [0.01, 1.0, 0.8, 1.4], ... % tissue OPs [ua, us', g, n] 
           [0.0, 1e-10, 1.0, 1.0] ... % air optical properties
       } ... 
   ); 
   

   The vector with comment "tissue OPs" specifies μ_a mm^-1, μ_s' mm^-1, g and n of the tissue slab. Fluence is calculated in a cylindrical coordinate system (ρ, z) for increasing number of photons launched from the source, N = 10, 100, 1000, 10000.

2. Configure an Analog simulation by setting "options.AbsorptionWeightingType" to 'Analog' (Line 17).

     options.AbsorptionWeightingType = 'Analog'; 
   

3. Set the tissue absorption to μ_a=0.01 mm^-1.
4. Execute Example 1a by right clicking on the cell so that it turns a pale yellow color and selecting "Evaluate Current Cell" or typing Ctrl+Enter.
5. How do the fluence results change as N is increased?
6. The relative error (RE) for the fluence is determined by dividing the standard deviation by the mean. For N=10,000, why does the relative error increase with distance from the source?
7. Change the Monte Carlo estimator to Discrete Absorption Weighting (DAW) by setting "options.AbsorptionWeightingType" to 'Discrete'.

     options.AbsorptionWeightingType = 'Discrete'; 
   

   Execute the Example 1a cell with this new setting.

8. How do the results using Analog and DAW compare for a given number of launched photons N?
9. For N=10000 which estimator provides fluence results with the least RE?
10. Now move down to Example 1b. This cell plots the difference in the relative error (RE) between Analog and DAW estimators.
11. Execute this cell. In what spatial regions does DAW show smaller RE than Analog? Can you explain why?

II. Spatially-Resolved Reflectance Predictions for Analog versus CAW Simulations

Goal: This exercise compares error estimates of spatially-resolved reflectance using Analog versus Continuous Absorption Weighting (CAW) simulations.

1. Example 2 sets up two Monte Carlo simulations using: 1) Analog, and 2) Continuous Absorption Weighting (CAW). The source is collimated and the tissue has optical properties [μ_a, μ'_s, g, n] = [0.01, 1.0, 0.8, 1.4]. Spatially-resolved reflectance is calculated in ρ bins. The number of photons launched for each simulation is N=10,000.
2. Set the tissue absorption to μ_a=0.01 mm^-1

      [0.01, 1.0, 0.8, 1.4], ... % tissue OPs [ua, us', g, n]
   

   and execute this cell.

3. What span of source-detector separations (ρ) does Analog have smaller standard deviations (SDs) than CAW? Where does CAW have smaller SDs? Can you explain why?
4. Note the time each took to execute in the Matlab command window. Why does Analog run so much faster?
5. Determine the Analog relative error in the last ρ bin (ρ=9.9 mm) by typing "d1SD(100)/d1.Mean(100)" in the matlab command window. Do the same for CAW "d2SD(100)/d2.Mean(100)". Which method, Analog or CAW, has the smaller relative error?
6. The efficiency of an estimator is calculated as Efficiency =1/(R^2 T) where R is the relative error and T is the time it took for the simulation to execute. Which estimator has the larger efficiency for detector estimates at ρ=9.9 mm?
7. Set the tissue absorption to μ_a=0.1 mm^-1

      [0.1, 1.0, 0.8, 1.4], ... % tissue OPs [ua, us', g, n]
   

   and execute the cell again.

8. Looking at the results for μ_a=0.01 mm^-1 with μ_a=0.1 mm^-1, how do the two methods RE compare as absorption is increased?