Skip to content

Commit 5258bf6

Browse files
lmalenfantlmalenfant
committed
Updated the Python code and added a Python file for the impulse beam demo
1 parent d7bd4cc commit 5258bf6

2 files changed

Lines changed: 114 additions & 4 deletions

File tree

‎monte-carlo/r-of-rho-impulse-beam.py‎

Lines changed: 111 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,111 @@
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()

‎monte-carlo/r-of-rho.py‎

Lines changed: 3 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,11 @@
11
# This is an example of python code using VTS to plot R(rho) using MCCL
22
#
33
# Import the Operating System so we can access the files for the VTS library
4-
import libraries
54
from pythonnet import load
65
load('coreclr')
76
import clr
87
import os
9-
file = './libraries/Vts.dll'
8+
file = '../libraries/Vts.dll'
109
print('Does this filepath exist?', os.path.isfile(file))
1110
clr.AddReference(os.path.abspath(file))
1211

@@ -41,8 +40,8 @@
4140
detectors[0] = detectorInput
4241

4342
simulationInput = SimulationInput()
44-
simulationInput.N=1000
45-
simulationInput.DetectorInputs= detectors
43+
simulationInput.N = 1000
44+
simulationInput.DetectorInputs = detectors
4645

4746
# create the simulation
4847
simulation = MonteCarloSimulation(simulationInput)

0 commit comments

Comments
 (0)