Skip to content

Commit 11c6a40

Browse files
lmalenfantlmalenfant
authored
Merge pull request #21 from VirtualPhotonics/feature/19-add-a-fluence-sample
2 parents 50d58ee + 9458d85 commit 11c6a40

3 files changed

Lines changed: 309 additions & 0 deletions

File tree

‎forward-solvers/fluence-of-rho-and-z-two-layer.py‎

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
# This is an example of python code using VTS to Compute fluence for a two-layer
2+
# medium as a function of radial extent and depth at a given set of optical properties.
3+
# A two-layer SDA forward solver is used to compute the fluence with optical properties
4+
# defined in opRegions[0] opRegions[1], and top layer thickness defined in topLayerThickness [mm].
5+
# The fluence as a function of rho and z is determined and when displayed, it is mirrored
6+
# to show full fluence.
7+
#
8+
# Import PythonNet
9+
from pythonnet import load
10+
load('coreclr')
11+
import clr
12+
# Import the Operating System so we can access the files for the VTS library
13+
import os
14+
file = '../libraries/Vts.dll'
15+
clr.AddReference(os.path.abspath(file))
16+
import numpy as np
17+
import plotly.graph_objects as go
18+
import plotly.express as px
19+
from Vts import *
20+
from Vts.Common import *
21+
from Vts.Extensions import *
22+
from Vts.Modeling.Optimizers import *
23+
from Vts.Modeling.ForwardSolvers import *
24+
from Vts.SpectralMapping import *
25+
from Vts.Factories import *
26+
from Vts.MonteCarlo import *
27+
from Vts.MonteCarlo.Sources import *
28+
from Vts.MonteCarlo.Tissues import *
29+
from Vts.MonteCarlo.Detectors import *
30+
from Vts.MonteCarlo.Factories import *
31+
from Vts.MonteCarlo.PhotonData import *
32+
from Vts.MonteCarlo.PostProcessing import *
33+
from System import Array, Double, Object, Func, Math
34+
clr.AddReference("System.Core")
35+
from System.Linq import Enumerable
36+
37+
solver = TwoLayerSDAForwardSolver()
38+
solver.SourceConfiguration = SourceConfiguration.Distributed
39+
topLayerThickness = 5
40+
opRegions = Array.CreateInstance(IOpticalPropertyRegion, 2)
41+
opRegions[0] = LayerOpticalPropertyRegion(DoubleRange(0, topLayerThickness, 2), OpticalProperties(0.1, 1, 0.8, 1.4))
42+
opRegions[1] = LayerOpticalPropertyRegion(DoubleRange(topLayerThickness, Double.PositiveInfinity, 2), OpticalProperties(0.01, 1, 0.8, 1.4))
43+
# Create the DoubleRange instance
44+
rhos_range = DoubleRange(0.1, 19.9, 100) # range of s-d separations in mm
45+
46+
# Convert to .NET array
47+
rho_delta = rhos_range.GetDelta()
48+
print(rho_delta)
49+
# Start at 0.1, increment by 0.2, 100 elements
50+
rhos = 0.1 + rho_delta * np.arange(100)
51+
print(rhos)
52+
53+
zs_range = DoubleRange(0.1, 19.9, 100) # range of depths in mm
54+
55+
# Convert to .NET array
56+
zs_delta = zs_range.GetDelta()
57+
print(zs_delta)
58+
# Start at 0.1, increment by delta, 100 elements
59+
zs = 0.1 + zs_delta * np.arange(100)
60+
print(zs)
61+
62+
arrayOP = Array[Array[IOpticalPropertyRegion]]([opRegions])
63+
# predict the tissue's fluence(rho, z) for the given optical properties
64+
fluenceOfRhoAndZ = solver.FluenceOfRhoAndZ(arrayOP, rhos, zs );
65+
66+
allRhos = np.concatenate((-rhos[::-1], rhos))
67+
68+
# log transform
69+
log_fluence = [Math.Log(f) for f in fluenceOfRhoAndZ]
70+
71+
size = len(zs)
72+
# split into rows
73+
fluenceRowsToPlot = np.array([log_fluence[i:i+size] for i in range(0, len(log_fluence), size)])
74+
75+
# reverse and concatenate
76+
allFluenceRowsToPlot = np.concatenate((fluenceRowsToPlot[::-1], fluenceRowsToPlot))
77+
78+
# Heatmap function to convert the data into a heat map
79+
def heatmap(values, x, y, x_label="", y_label="", title=""):
80+
"""Create a heatmap chart."""
81+
# values should be a 2D array-like (list of lists or 2D numpy array)
82+
fig = go.Figure(data=go.Heatmap(
83+
z=values,
84+
x=x,
85+
y=y,
86+
transpose=True,
87+
colorscale='Hot',
88+
colorbar=dict(title=title)
89+
))
90+
fig.update_layout(
91+
title=title,
92+
xaxis_title=x_label,
93+
yaxis_title=y_label,
94+
yaxis_autorange='reversed'
95+
)
96+
return fig
97+
98+
fluenceChart = heatmap(allFluenceRowsToPlot.tolist(), allRhos.tolist(), list(zs), "ρ [mm]", "z [mm]", "log(Φ(ρ, z) [mm-2])")
99+
fluenceChart.add_hline(y=topLayerThickness, line_dash="dash", line_color="white", line_width=2)
100+
fluenceChart.show(renderer="browser")

‎forward-solvers/phd-of-rho-and-z-compute-fluence.py‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
# This is an example of python code using VTS to Compute photon hitting density for homogeneous
2+
# medium at a given set of optical properties opRegions[0].
3+
# This sample uses ComputeFluence in place of calling FluenceOfRhoAndZ on the forward solver object
4+
# and it uses a distributed point source SDA Forward Solver
5+
#
6+
# Import PythonNet
7+
from pythonnet import load
8+
load('coreclr')
9+
import clr
10+
# Import the Operating System so we can access the files for the VTS library
11+
import os
12+
file = '../libraries/Vts.dll'
13+
clr.AddReference(os.path.abspath(file))
14+
import numpy as np
15+
import plotly.graph_objects as go
16+
import plotly.express as px
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+
from System import Array, Double, Object, Func, Math
32+
clr.AddReference("System.Core")
33+
from System.Linq import Enumerable
34+
35+
solver = DistributedPointSourceSDAForwardSolver()
36+
37+
topLayerThickness = 5
38+
opRegions = Array.CreateInstance(IOpticalPropertyRegion, 1)
39+
opRegions[0] = LayerOpticalPropertyRegion(DoubleRange(0, topLayerThickness, 2), OpticalProperties(0.1, 1, 0.8, 1.4))
40+
# Create the DoubleRange instance
41+
rhos_range = DoubleRange(0.1, 19.9, 100) # range of s-d separations in mm
42+
43+
# Convert to .NET array
44+
rho_delta = rhos_range.GetDelta()
45+
print(rho_delta)
46+
# Start at 0.1, increment by 0.2, 100 elements
47+
rhos = 0.1 + rho_delta * np.arange(100)
48+
print(rhos)
49+
50+
zs_range = DoubleRange(0.1, 19.9, 100) # range of depths in mm
51+
52+
# Convert to .NET array
53+
zs_delta = zs_range.GetDelta()
54+
print(zs_delta)
55+
# Start at 0.1, increment by delta, 100 elements
56+
zs = 0.1 + zs_delta * np.arange(100)
57+
print(zs)
58+
59+
allRhos = np.concatenate((-rhos[::-1], rhos))
60+
61+
opRegionsArray = Array[Array[IOpticalPropertyRegion]]([opRegions])
62+
# predict the tissue's fluence(rho, z) for the given optical properties
63+
independentAxes = Array.CreateInstance(IndependentVariableAxis, 1)
64+
independentAxes[0] = IndependentVariableAxis.Z
65+
independentValues = Array.CreateInstance(Array[Double], 2)
66+
independentValues[0] = Array[Double](allRhos.tolist())
67+
independentValues[1] = Array[Double](zs.tolist())
68+
# Call the static method ComputeFluence in ComputationFactory to get the fluence data
69+
fluenceOfRhoAndZ = ComputationFactory.ComputeFluence(solver, FluenceSolutionDomainType.FluenceOfRhoAndZ, independentAxes, independentValues, opRegions, Array[Double](allRhos.tolist()))
70+
71+
#PHD
72+
sourceDetectorSeparation = 10
73+
opArray = Array.CreateInstance(OpticalProperties, 2)
74+
opArray[0] = OpticalProperties(0.1, 1, 0.8, 1.4)
75+
opArray[1] = OpticalProperties(0.01, 1, 0.8, 1.4)
76+
77+
phdOfRhoAndZ = ComputationFactory.GetPHD(ForwardSolverType.TwoLayerSDA, fluenceOfRhoAndZ, sourceDetectorSeparation, opArray, Array[Double](allRhos.tolist()), Array[Double](zs.tolist()))
78+
79+
# log transform
80+
log_phd = [Math.Log(f) for f in phdOfRhoAndZ]
81+
82+
size = len(zs)
83+
# split into rows
84+
phdRowsToPlot = np.array([log_phd[i:i+size] for i in range(0, len(log_phd), size)])
85+
86+
# Heatmap function to convert the data into a heat map
87+
def heatmap(values, x, y, x_label="", y_label="", title=""):
88+
"""Create a heatmap chart."""
89+
# values should be a 2D array-like (list of lists or 2D numpy array)
90+
fig = go.Figure(data=go.Heatmap(
91+
z=values,
92+
x=x,
93+
y=y,
94+
transpose=True,
95+
colorscale='Hot',
96+
colorbar=dict(title=title)
97+
))
98+
fig.update_layout(
99+
title=title,
100+
xaxis_title=x_label,
101+
yaxis_title=y_label,
102+
yaxis_autorange='reversed'
103+
)
104+
return fig
105+
106+
fluenceChart = heatmap(phdRowsToPlot.tolist(), allRhos.tolist(), list(zs), "ρ [mm]", "z [mm]", "log(phd(ρ, z) [mm-2])")
107+
fluenceChart.show(renderer="browser")

‎forward-solvers/phd-of-rho-and-z-two-layer.py‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
# This is an example of python code using VTS to Compute photon hitting density for a two-layer
2+
# medium at a given set of optical properties. A two-layer SDA forward solver is used to
3+
# compute the photon hitting density with optical properties defined in opRegions[0] and
4+
# opRegions[1], and top layer thickness defined in topLayerThickness [mm].
5+
#
6+
# Import PythonNet
7+
from pythonnet import load
8+
load('coreclr')
9+
import clr
10+
# Import the Operating System so we can access the files for the VTS library
11+
import os
12+
file = '../libraries/Vts.dll'
13+
clr.AddReference(os.path.abspath(file))
14+
import numpy as np
15+
import plotly.graph_objects as go
16+
import plotly.express as px
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+
from System import Array, Double, Object, Func, Math
32+
clr.AddReference("System.Core")
33+
from System.Linq import Enumerable
34+
35+
solver = TwoLayerSDAForwardSolver()
36+
solver.SourceConfiguration = SourceConfiguration.Distributed
37+
topLayerThickness = 5
38+
opRegions = Array.CreateInstance(IOpticalPropertyRegion, 2)
39+
opRegions[0] = LayerOpticalPropertyRegion(DoubleRange(0, topLayerThickness, 2), OpticalProperties(0.1, 1, 0.8, 1.4))
40+
opRegions[1] = LayerOpticalPropertyRegion(DoubleRange(topLayerThickness, Double.PositiveInfinity, 2), OpticalProperties(0.01, 1, 0.8, 1.4))
41+
# Create the DoubleRange instance
42+
rhos_range = DoubleRange(0.1, 19.9, 100) # range of s-d separations in mm
43+
44+
# Convert to .NET array
45+
rho_delta = rhos_range.GetDelta()
46+
print(rho_delta)
47+
# Start at 0.1, increment by 0.2, 100 elements
48+
rhos = 0.1 + rho_delta * np.arange(100)
49+
print(rhos)
50+
51+
zs_range = DoubleRange(0.1, 19.9, 100) # range of depths in mm
52+
53+
# Convert to .NET array
54+
zs_delta = zs_range.GetDelta()
55+
print(zs_delta)
56+
# Start at 0.1, increment by delta, 100 elements
57+
zs = 0.1 + zs_delta * np.arange(100)
58+
print(zs)
59+
60+
allRhos = np.concatenate((-rhos[::-1], rhos))
61+
opRegionsArray = Array[Array[IOpticalPropertyRegion]]([opRegions])
62+
# predict the tissue's fluence(rho, z) for the given optical properties
63+
fluenceOfRhoAndZ = solver.FluenceOfRhoAndZ(opRegionsArray, allRhos, zs );
64+
65+
#PHD
66+
sourceDetectorSeparation = 10
67+
opArray = Array.CreateInstance(OpticalProperties, 2)
68+
opArray[0] = OpticalProperties(0.1, 1, 0.8, 1.4)
69+
opArray[1] = OpticalProperties(0.01, 1, 0.8, 1.4)
70+
71+
phdOfRhoAndZ = ComputationFactory.GetPHD(solver, fluenceOfRhoAndZ, sourceDetectorSeparation, opArray, Array[Double](allRhos.tolist()), Array[Double](zs.tolist()))
72+
73+
# log transform
74+
log_fluence = [Math.Log(f) for f in phdOfRhoAndZ]
75+
76+
size = len(zs)
77+
# split into rows
78+
fluenceRowsToPlot = np.array([log_fluence[i:i+size] for i in range(0, len(log_fluence), size)])
79+
80+
# Heatmap function to convert the data into a heat map
81+
def heatmap(values, x, y, x_label="", y_label="", title=""):
82+
"""Create a heatmap chart."""
83+
# values should be a 2D array-like (list of lists or 2D numpy array)
84+
fig = go.Figure(data=go.Heatmap(
85+
z=values,
86+
x=x,
87+
y=y,
88+
transpose=True,
89+
colorscale='Hot',
90+
colorbar=dict(title=title)
91+
))
92+
fig.update_layout(
93+
title=title,
94+
xaxis_title=x_label,
95+
yaxis_title=y_label,
96+
yaxis_autorange='reversed'
97+
)
98+
return fig
99+
100+
fluenceChart = heatmap(fluenceRowsToPlot.tolist(), allRhos.tolist(), list(zs), "ρ [mm]", "z [mm]", "log(phd(ρ, z) [mm-2])")
101+
fluenceChart.add_hline(y=topLayerThickness, line_dash="dash", line_color="white", line_width=2)
102+
fluenceChart.show(renderer="browser")

0 commit comments

Comments
 (0)