1+ # This is an example of python code using VTS to plot R(rho) using MCCL
2+ #
3+ # Import the Operating System so we can access the files for the VTS library
4+ from pythonnet import load
5+ load ('coreclr' )
6+ import clr
7+ import os
8+ file = '../libraries/Vts.dll'
9+ print ('Does this filepath exist?' , os .path .isfile (file ))
10+ clr .AddReference (os .path .abspath (file ))
11+
12+ print ('Import numpy' )
13+ import numpy as np
14+ print ('Import pyplot' )
15+ import matplotlib .pyplot as plt
16+ print ('Import Vts' )
17+ from Vts import *
18+ from Vts .Common import *
19+ from Vts .Extensions import *
20+ from Vts .Modeling .Optimizers import *
21+ from Vts .Modeling .ForwardSolvers import *
22+ from Vts .SpectralMapping import *
23+ from Vts .Factories import *
24+ from Vts .MonteCarlo import *
25+ from Vts .MonteCarlo .Sources import *
26+ from Vts .MonteCarlo .Tissues import *
27+ from Vts .MonteCarlo .Detectors import *
28+ from Vts .MonteCarlo .Factories import *
29+ from Vts .MonteCarlo .PhotonData import *
30+ from Vts .MonteCarlo .PostProcessing import *
31+ print ('Import System' )
32+ from System import Array , Double
33+
34+ # SimulationInput defines the simulation. I think the default is collimated point source illumination normal to the surface.
35+ simulationInput = SimulationInput ()
36+ simulationInput .N = 10000
37+
38+ # We add four detectors. One to monitor the reflected excitance as a function of radius, and three to monitor the overall reflection and transmission.
39+ # add detectors to Simulation object
40+ detectors = Array .CreateInstance (IDetectorInput , 4 )
41+
42+ detectors [0 ] = ROfRhoDetectorInput ()
43+ detectors [0 ].Rho = DoubleRange (start = 0 , stop = 5 , number = 401 )
44+ detectors [1 ] = RDiffuseDetectorInput ()
45+ detectors [2 ] = RSpecularDetectorInput ()
46+ detectors [3 ] = TDiffuseDetectorInput ()
47+
48+ simulationInput .DetectorInputs = detectors
49+
50+ # This is an index matched sample that has scattering of 0.99/mm and absorption of 0.01/mm. It is 10mm thick.
51+ d = 10 # mm
52+ mua = 0.01 # mm⁻¹
53+ musp = 0.99 # mm⁻¹
54+ g = 0 # scattering anisotropy
55+ n = 1 # matched index of refraction
56+
57+ regions = Array .CreateInstance (ITissueRegion , 3 )
58+
59+ regions [0 ] = LayerTissueRegion (zRange = DoubleRange (Double .NegativeInfinity , 0.0 ), op = OpticalProperties (mua = 0.0 , musp = 1E-10 , g = 1.0 , n = 1.0 ))
60+ regions [1 ] = LayerTissueRegion (zRange = DoubleRange (0.0 , d ), op = OpticalProperties (mua = mua , musp = musp , g = g , n = n ))
61+ regions [2 ] = LayerTissueRegion (zRange = DoubleRange (d , Double .PositiveInfinity ), op = OpticalProperties (0.0 , 1E-10 , 1.0 , 1.0 ))
62+
63+ simulationInput .TissueInput = MultiLayerTissueInput (regions )
64+
65+ # create the simulation
66+ simulation = MonteCarloSimulation (simulationInput )
67+
68+ # run the simulation
69+ simulationOutput = simulation .Run ()
70+
71+ # extract the data
72+ rDiffuse = Array .CreateInstance (RDiffuseDetector , 1 )
73+ rDiffuse [0 ] = simulationOutput .ResultsDictionary ["RDiffuse" ]
74+
75+ tDiffuse = Array .CreateInstance (TDiffuseDetector , 1 )
76+ tDiffuse [0 ] = simulationOutput .ResultsDictionary ["TDiffuse" ]
77+
78+ rSpecular = Array .CreateInstance (RSpecularDetector , 1 )
79+ rSpecular [0 ] = simulationOutput .ResultsDictionary ["RSpecular" ]
80+
81+ print ("specular R = %.5f" % rSpecular [0 ].Mean )
82+ print (" diffuse R = %.5f" % rDiffuse [0 ].Mean )
83+ print (" diffuse T = %.5f" % tDiffuse [0 ].Mean )
84+
85+ # plot the results using Pyplot
86+ detectorResults = Array .CreateInstance (ROfRhoDetector , 1 )
87+ detectorResults [0 ] = simulationOutput .ResultsDictionary ["ROfRho" ]
88+
89+ reflectance = np .array ([r for r in detectorResults [0 ].Mean ])
90+ edges = np .array ([mp for mp in detectorResults [0 ].Rho ])[:- 1 ]
91+
92+ plt .figure (figsize = (8 ,4.5 ))
93+ plt .plot (edges , reflectance , 'ob' , markersize = 2 )
94+ plt .plot (- edges , reflectance , 'ob' , markersize = 2 )
95+ plt .xlabel ('ρ [mm]' )
96+ plt .ylabel ('R(ρ) [W/mm²]' )
97+ plt .title ('1 W Impulse Beam' )
98+
99+ text_args = {
100+ 'ha' : 'left' , # Horizontal alignment
101+ 'va' : 'top' , # Vertical alignment
102+ 'transform' : plt .gca ().transAxes , # Coordinate system transformation
103+ 'fontsize' : 9 # Font size
104+ }
105+
106+ plt .text (0.8 , 0.95 , 'n=%.3f' % n , ** text_args )
107+ plt .text (0.8 , 0.90 , r'$\mu_a$=%.2f mm⁻¹' % mua , ** text_args )
108+ plt .text (0.8 , 0.85 , r"${\mu_s}'$=%.2f mm⁻¹" % musp , ** text_args )
109+ plt .text (0.8 , 0.80 , 'g=%.2f' % g , ** text_args )
110+ plt .text (0.8 , 0.75 , 'd=%.1f mm' % d , ** text_args )
111+ plt .show ()
0 commit comments