Skip to content

Commit 50d58ee

Browse files
hayakawa16hayakawa16
authored
Merge pull request #20 from VirtualPhotonics/feature/18-add-inverse-solution-scripts
Feature/18 add inverse solution scripts
2 parents 84b87a0 + 7bd1f25 commit 50d58ee

3 files changed

Lines changed: 446 additions & 0 deletions

File tree

‎inverse-solutions/r-of-fx-multi-op-prop-inversion.py‎

Lines changed: 176 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,176 @@
1+
# This is an example of python code using VTS to provide inverse solution
2+
# for R(fx) to find chromophore concentrations [Hb HbO2] and power law
3+
# coefficients [A b] using fixed fx=[0 0.2]/mm and wavelengths=[650:50:1000]nm.
4+
# Scaled Monte Carlo with Nurbs forward solver provides the simulated
5+
# measured data and PointSourceSDA provides the model used during the inversion.
6+
# The optimization is performed by a python library scipy.
7+
#
8+
# Import the Operating System so we can access the files for the VTS library
9+
from pythonnet import load
10+
load('coreclr')
11+
import clr
12+
import os
13+
file = '../libraries/Vts.dll'
14+
clr.AddReference(os.path.abspath(file))
15+
import numpy as np
16+
import plotly.graph_objects as go
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, Object
32+
# Setup wavelengths in visible and NIR spectral regimes
33+
wavelengths = Array.CreateInstance(float, 8)
34+
for i in range(0, len(wavelengths)):
35+
wavelengths[i] = 650.0 + 50 * i
36+
# Setup fxs
37+
fxs = [0.0, 0.2]
38+
# Define parameters to fit [Hb,HbO2,A,b]
39+
measuredData = [28.4, 22.4, 1.2, 1.42]
40+
# Construct a scatter
41+
scattererMeasuredData = PowerLawScatterer(measuredData[2], measuredData[3])
42+
# Setup the values for the measured data
43+
chromophoresMeasuredData = Array.CreateInstance(IChromophoreAbsorber, 2)
44+
chromophoresMeasuredData[0] = ChromophoreAbsorber(ChromophoreType.HbO2, measuredData[0])
45+
chromophoresMeasuredData[1] = ChromophoreAbsorber(ChromophoreType.Hb, measuredData[1])
46+
# len(opsMeasured)=8 per number of wavelengths
47+
opsMeasured = Tissue(chromophoresMeasuredData, scattererMeasuredData, "", n=1.4).GetOpticalProperties(wavelengths)
48+
# Create measurements using Nurbs-based white Monte Carlo forward solver
49+
measurementForwardSolver = NurbsForwardSolver()
50+
# len(measuredROfFx)=(#fxs)x(#wavelengths) flattened
51+
rOfFxMeasured=np.concatenate(
52+
[np.array(measurementForwardSolver.ROfFx(opsMeasured, fxs[0]), dtype=float),
53+
np.array(measurementForwardSolver.ROfFx(opsMeasured, fxs[1]), dtype=float)])
54+
# Create a forward solver as a model function for inversion
55+
forwardSolverForInversion = PointSourceSDAForwardSolver()
56+
#forwardSolverForInversion = NurbsForwardSolver() # results improved but inverse crime!
57+
58+
# Declare local forward reflectance function that computes reflectance
59+
# from chromophores and scatterer values
60+
# valuesSought = [Hb, HbO2, A, b]
61+
def CalculateReflectanceVsWavelengthFromChromophoreConcAndScatterer(
62+
valuesSought, wavelengths, fxs, forwardSolver):
63+
# Create a forward solve model function to solve inverse
64+
forwardSolverForInversion = PointSourceSDAForwardSolver()
65+
# Create an array of chromophore absorbers based on values
66+
chromophoresLocal = Array.CreateInstance(IChromophoreAbsorber, 2)
67+
chromophoresLocal[0] = ChromophoreAbsorber(ChromophoreType.HbO2, valuesSought[0])
68+
chromophoresLocal[1] = ChromophoreAbsorber(ChromophoreType.Hb, valuesSought[1])
69+
# Create a scatterr based on values
70+
scattererLocal = PowerLawScatterer(valuesSought[2], valuesSought[3])
71+
# Compose local tissue to obtain optical properties
72+
opsLocal = Tissue(chromophoresLocal, scattererLocal, "", n=1.4).GetOpticalProperties(wavelengths)
73+
print("iter:[Hb HbO2 A b]=[%5.3f %5.3f %5.3f %5.3f]" % (
74+
valuesSought[0], valuesSought[1], valuesSought[2],valuesSought[3]))
75+
# Compute reflectance for local absorbers
76+
modelDataLocal=np.concatenate(
77+
[np.array(forwardSolver.ROfFx(opsLocal, fxs[0]), dtype=float),
78+
np.array(forwardSolver.ROfFx(opsLocal, fxs[1]), dtype=float)])
79+
return modelDataLocal
80+
81+
# func for residual
82+
def residual(valuesSought, wavelengths, fxs, measuredROfFx, forwardSolver):
83+
predictedROfFx= CalculateReflectanceVsWavelengthFromChromophoreConcAndScatterer(
84+
valuesSought, wavelengths, fxs, forwardSolver)
85+
difference = Array.CreateInstance(float,len(wavelengths)*len(fxs))
86+
for i in range(0, len(fxs)*len(wavelengths)-1):
87+
difference[i] = predictedROfFx[i] - measuredROfFx[i]
88+
return difference
89+
90+
# Run the inversion: set up initial guess
91+
initialGuess = [18.0, 30.0, 0.8, 1.6]
92+
chromophoresInitialGuess = Array.CreateInstance(IChromophoreAbsorber, 2)
93+
chromophoresInitialGuess[0] = ChromophoreAbsorber(ChromophoreType.HbO2, initialGuess[0])
94+
chromophoresInitialGuess[1] = ChromophoreAbsorber(ChromophoreType.Hb, initialGuess[1])
95+
scattererInitialGuess = PowerLawScatterer(initialGuess[2], initialGuess[3])
96+
# Compose tissue for initial guess data to obtain OPs
97+
opsInitialGuess = Tissue(chromophoresInitialGuess, scattererInitialGuess, "", n=1.4).GetOpticalProperties(wavelengths)
98+
rOfFxInitialGuess=np.concatenate(
99+
[np.array(forwardSolverForInversion.ROfFx(opsInitialGuess, fxs[0]), dtype=float),
100+
np.array(forwardSolverForInversion.ROfFx(opsInitialGuess, fxs[1]), dtype=float)])
101+
# Run the levenberg-marquardt inversion
102+
from scipy.optimize import least_squares
103+
fit = least_squares(
104+
residual,
105+
initialGuess,
106+
ftol=1e-9, xtol=1e-9, max_nfev=10000, # max_nfev needs to be integer
107+
args=(wavelengths, fxs, rOfFxMeasured, forwardSolverForInversion),
108+
method='lm')
109+
# Calculate final reflectance from model at fit values
110+
chromophoresFit = Array.CreateInstance(IChromophoreAbsorber, 2)
111+
chromophoresFit[0] = ChromophoreAbsorber(ChromophoreType.HbO2, fit.x[0])
112+
chromophoresFit[1] = ChromophoreAbsorber(ChromophoreType.Hb, fit.x[1])
113+
scattererFit = PowerLawScatterer(fit.x[2], fit.x[3])
114+
# Compose tissue for fit data to obtain OPs
115+
opsFit= Tissue(chromophoresFit, scattererFit, "", n=1.4).GetOpticalProperties(wavelengths)
116+
rOfFxFit=np.concatenate(
117+
[np.array(forwardSolverForInversion.ROfFx(opsFit, fxs[0]), dtype=float),
118+
np.array(forwardSolverForInversion.ROfFx(opsFit, fxs[1]), dtype=float)])
119+
# plot Reflectance: flattened so have to separate
120+
chart1 = go.Figure()
121+
xLabel = "wavelength [nm]"
122+
yLabel = "R(wavelength)"
123+
wvs = [w for w in wavelengths]
124+
# plot measured data first fx first
125+
measR= [m for m in rOfFxMeasured]
126+
midpoint=len(measR) // 2
127+
chart1.add_trace(go.Scatter(x=wvs, y=measR[:midpoint], mode='markers', name='measured data: fx1'))
128+
chart1.add_trace(go.Scatter(x=wvs, y=measR[midpoint:], mode='markers', name='measured data: fx2'))
129+
# plot initial guess data
130+
igR = [i for i in rOfFxInitialGuess]
131+
chart1.add_trace(go.Scatter(x=wvs, y=igR[:midpoint], mode='markers', name='initial guess: fx1'))
132+
chart1.add_trace(go.Scatter(x=wvs, y=igR[midpoint:], mode='markers', name='initial guess: fx2'))
133+
# plot fit: need to organize by fx
134+
convR = [f for f in rOfFxFit]
135+
chart1.add_trace(go.Scatter(x=wvs, y=convR[:midpoint], mode='lines', name='converged: fx1'))
136+
chart1.add_trace(go.Scatter(x=wvs, y=convR[midpoint:], mode='lines', name='converged: fx2'))
137+
chart1.update_layout( title="ROfFx (inverse solution for chromophore concentrations, multiple wavelengths, multiple fx)", xaxis_title=xLabel, yaxis_title=yLabel)
138+
chart1.show(renderer="browser")
139+
# plot Mus': flattened so have to separate
140+
chart2 = go.Figure()
141+
xLabel = "wavelength [nm]"
142+
yLabel = "us'(wavelength)"
143+
wvs = [w for w in wavelengths]
144+
# plot measured data
145+
scattererMeasuredDataMusp= np.zeros(len(wavelengths),dtype=float)
146+
for i in range(0, len(wavelengths)):
147+
scattererMeasuredDataMusp[i]=opsMeasured[i].Musp
148+
measMusp = [m for m in scattererMeasuredDataMusp]
149+
chart2.add_trace(go.Scatter(x=wvs, y=measMusp, mode='markers', name='measured data'))
150+
# plot initial guess data
151+
scattererInitialGuessMusp= np.zeros(len(wavelengths),dtype=float)
152+
for i in range(0, len(wavelengths)):
153+
scattererInitialGuessMusp[i]=opsInitialGuess[i].Musp
154+
igMusp = [i for i in scattererInitialGuessMusp]
155+
chart2.add_trace(go.Scatter(x=wvs, y=igMusp, mode='markers', name='initial guess'))
156+
# plot fit
157+
scattererFitMusp= np.zeros(len(wavelengths),dtype=float)
158+
for i in range(0, len(wavelengths)):
159+
scattererFitMusp[i]=opsFit[i].Musp
160+
convMusp = [f for f in scattererFitMusp]
161+
chart2.add_trace(go.Scatter(x=wvs, y=convMusp, mode='lines', name='converged'))
162+
chart2.update_layout( title="ROfFx (inverse solution for chromophore concentrations, multiple wavelengths, multiple fx)", xaxis_title=xLabel, yaxis_title=yLabel)
163+
chart2.show(renderer="browser")
164+
# output results
165+
print("Meas = [%5.3f %5.3f %5.3f %5.3f]" % (
166+
measuredData[0], measuredData[1], measuredData[2], measuredData[3]))
167+
print("IG = [%5.3f %5.3f %5.3f %5.3f] Chi2=%5.3e" % (
168+
initialGuess[0], initialGuess[1], initialGuess[2], initialGuess[3],
169+
np.dot(rOfFxMeasured-rOfFxInitialGuess,rOfFxMeasured-rOfFxInitialGuess)))
170+
print("Conv = [%5.3f %5.3f %5.3f %5.3f] Chi2=%5.3e" % (fit.x[0], fit.x[1], fit.x[2], fit.x[3],
171+
np.dot(rOfFxMeasured-rOfFxFit,rOfFxMeasured-rOfFxFit)))
172+
print("error = [%3.2f %3.2f %3.2f %3.2f]%%" % (
173+
100*abs((measuredData[0]-fit.x[0])/measuredData[0]),
174+
100*abs((measuredData[1]-fit.x[1])/measuredData[1]),
175+
100*abs((measuredData[2]-fit.x[2])/measuredData[2]),
176+
100*abs((measuredData[3]-fit.x[3])/measuredData[3])))

‎inverse-solutions/r-of-rho-multi-op-prop-inversion-using-csharp-solve.py‎

Lines changed: 135 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,135 @@
1+
# This is an example of python code using VTS to provide inverse solution
2+
# for R(rho) to find chromophore concentrations of [Hb HbO2 H2O], with
3+
# fixed single rho=1mm, 13 wavelengths [400:50:1000]nm and scatterer
4+
# Power Law coefficients A=1.2, b=1.42.
5+
# Scaled Monte Carlo with Nurbs forward solver provides the simulated
6+
# measured data and PointSourceSDA provides the model used during the
7+
# inversion.
8+
# The optimization is performed by the Vts library MPFitLevenbergMarquardt
9+
# method Solve.
10+
#
11+
# Import the Operating System so we can access the files for the VTS library
12+
from pythonnet import load
13+
load('coreclr')
14+
import clr
15+
import os
16+
file = '../libraries/Vts.dll'
17+
clr.AddReference(os.path.abspath(file))
18+
import numpy as np
19+
import plotly.graph_objects as go
20+
from Vts import *
21+
from Vts.Common import *
22+
from Vts.Extensions import *
23+
from Vts.Modeling.Optimizers import *
24+
from Vts.Modeling.ForwardSolvers import *
25+
from Vts.SpectralMapping import *
26+
from Vts.Factories import *
27+
from Vts.MonteCarlo import *
28+
from Vts.MonteCarlo.Sources import *
29+
from Vts.MonteCarlo.Tissues import *
30+
from Vts.MonteCarlo.Detectors import *
31+
from Vts.MonteCarlo.Factories import *
32+
from Vts.MonteCarlo.PhotonData import *
33+
from Vts.MonteCarlo.PostProcessing import *
34+
from System import Array, Object, Func
35+
# Construct a scatterer
36+
scatterer = PowerLawScatterer(1.2, 1.42)
37+
# Setup wavelengths in visible and NIR spectral regimes
38+
wavelengths = Array.CreateInstance(float, 13)
39+
for i in range(0, len(wavelengths)):
40+
wavelengths[i] = 400.0 + 50 * i
41+
# Define rho
42+
rho = 1.0
43+
# Setup the values for the measured data
44+
measuredData = [70.0, 30.0, 0.8]
45+
chromophoresMeasuredData = Array.CreateInstance(IChromophoreAbsorber, 3)
46+
chromophoresMeasuredData[0] = ChromophoreAbsorber(ChromophoreType.HbO2, measuredData[0])
47+
chromophoresMeasuredData[1] = ChromophoreAbsorber(ChromophoreType.Hb, measuredData[1])
48+
chromophoresMeasuredData[2] = ChromophoreAbsorber(ChromophoreType.H2O, measuredData[2])
49+
opsMeasured = Tissue(chromophoresMeasuredData, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
50+
# Create measurements using Nurbs-based white Monte Carlo forward solver
51+
measurementForwardSolver = NurbsForwardSolver()
52+
rOfRhoMeasured = measurementForwardSolver.ROfRho(opsMeasured, rho)
53+
# Create a forward solver as a model function for inversion
54+
forwardSolverForInversion = PointSourceSDAForwardSolver()
55+
56+
# Declare local forward reflectance function that computes reflectance from chromophores
57+
# params=[wavelengths, rho, scatterer]
58+
def CalculateReflectanceFuncVsWavelengthFromChromophoreConcentration(
59+
chromophoreConcentration, params):
60+
# Create an array of chromophore absorbers based on values
61+
chromophoresLocal = Array.CreateInstance(IChromophoreAbsorber, 3)
62+
chromophoresLocal[0] = ChromophoreAbsorber(ChromophoreType.HbO2, chromophoreConcentration[0])
63+
chromophoresLocal[1] = ChromophoreAbsorber(ChromophoreType.Hb, chromophoreConcentration[1])
64+
chromophoresLocal[2] = ChromophoreAbsorber(ChromophoreType.H2O, chromophoreConcentration[2])
65+
# Compose local tissue to obtain optical properties
66+
opsLocal = Tissue(chromophoresLocal, params[2], "", n=1.4).GetOpticalProperties(params[0])
67+
print("iter:[HbO2,Hb,H2O]=[%3.2f %3.2f %3.2f]" % (
68+
chromophoreConcentration[0],chromophoreConcentration[1],chromophoreConcentration[2]))
69+
# Compute reflectance for local absorbers
70+
modelDataLocal = forwardSolverForInversion.ROfRho(opsLocal, params[1])
71+
return modelDataLocal
72+
73+
# Convert the Python function to a .NET Func delegate
74+
forward_func = Func[Array[float], Array[Object], Array[float]](CalculateReflectanceFuncVsWavelengthFromChromophoreConcentration)
75+
76+
# Run the inversion: set up initial guess
77+
initialGuess = [70.0, 30.0, 0.8]
78+
chromophoresInitialGuess = Array.CreateInstance(IChromophoreAbsorber, 3)
79+
chromophoresInitialGuess[0] = ChromophoreAbsorber(ChromophoreType.HbO2, initialGuess[0])
80+
chromophoresInitialGuess[1] = ChromophoreAbsorber(ChromophoreType.Hb, initialGuess[1])
81+
chromophoresInitialGuess[2] = ChromophoreAbsorber(ChromophoreType.H2O, initialGuess[2])
82+
# Compose tissue for initial guess data to obtain OPs
83+
opsInitialGuess = Tissue(chromophoresInitialGuess, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
84+
rOfRhoInitialGuess = forwardSolverForInversion.ROfRho(opsInitialGuess, rho)
85+
# Run the levenberg-marquardt inversion
86+
optimizer = MPFitLevenbergMarquardtOptimizer()
87+
initialGuessCopy = Array.CreateInstance(float, 3)
88+
initialGuessCopy = initialGuess
89+
parametersToFit = Array.CreateInstance(bool, 3)
90+
parametersToFit = [True, True, True]
91+
measuredDataWeight = Array.CreateInstance(float, len(rOfRhoMeasured))
92+
measuredDataWeight = [1] * len(measuredDataWeight)
93+
params = Array.CreateInstance(Object, 3)
94+
params = [wavelengths, rho, scatterer]
95+
96+
# try calling our LM
97+
fit = optimizer.Solve(initialGuessCopy, parametersToFit, rOfRhoMeasured, measuredDataWeight,
98+
forward_func, params)
99+
# Calculate final reflectance from model at fit values
100+
chromophoresFit = Array.CreateInstance(IChromophoreAbsorber, 3)
101+
chromophoresFit[0] = ChromophoreAbsorber(ChromophoreType.HbO2, fit[0])
102+
chromophoresFit[1] = ChromophoreAbsorber(ChromophoreType.Hb, fit[1])
103+
chromophoresFit[2] = ChromophoreAbsorber(ChromophoreType.H2O, fit[2])
104+
opsFit = Tissue(chromophoresFit, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
105+
rOfRhoFit= forwardSolverForInversion.ROfRho(opsFit, rho)
106+
# plot the results using Plotly
107+
xLabel = "wavelength [nm]"
108+
yLabel = "R(wavelength) [mm-2]"
109+
wvs = [w for w in wavelengths]
110+
# plot measured data
111+
meas = [m for m in rOfRhoMeasured]
112+
chart = go.Figure()
113+
chart.add_trace(go.Scatter(x=wvs, y=meas, mode='markers', name='measured data'))
114+
# plot initial guess data
115+
ig = [i for i in rOfRhoInitialGuess]
116+
chart.add_trace(go.Scatter(x=wvs, y=ig, mode='markers', name='initial guess'))
117+
# plot fit
118+
conv = [f for f in rOfRhoFit]
119+
chart.add_trace(go.Scatter(x=wvs, y=conv, mode='lines', name='converged'))
120+
chart.update_layout( title="ROfRho (inverse solution for chromophore concentrations, multiple wavelengths, single rho)", xaxis_title=xLabel, yaxis_title=yLabel)
121+
chart.show(renderer="browser")
122+
# output results
123+
print("Meas = [%5.3f %5.3f %5.3f]" % (measuredData[0], measuredData[1], measuredData[2]))
124+
print("IG = [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (
125+
initialGuess[0], initialGuess[1], initialGuess[2],
126+
np.dot(np.subtract(rOfRhoMeasured,rOfRhoInitialGuess),
127+
np.subtract(rOfRhoMeasured,rOfRhoInitialGuess))))
128+
print("Conv = [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (fit[0], fit[1], fit[2],
129+
np.dot(np.subtract(rOfRhoMeasured,rOfRhoFit),
130+
np.subtract(rOfRhoMeasured,rOfRhoFit))))
131+
print("error = [%5.3f %5.3f %5.3f]%%" % (
132+
100*abs((measuredData[0]-fit[0])/measuredData[0]),
133+
100*abs((measuredData[1]-fit[1])/measuredData[1]),
134+
100*abs((measuredData[2]-fit[2])/measuredData[2])))
135+

0 commit comments

Comments
 (0)