Skip to content

Commit b0172cc

Browse files
hayakawa16hayakawa16
committed
First draft of R(fx) inversion across 2 fxs and 8 wavelengths. Not running yet but didn't want to loose my work.
1 parent e6e571b commit b0172cc

1 file changed

Lines changed: 137 additions & 0 deletions

File tree

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

Lines changed: 137 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
# This is an example of python code using VTS to provide inverse solution
2+
# for R(fx) to find chromophore concentrations and power law coefficients.
3+
# The optimization uses the Vts library MPFitLevenbergMarquardt method Solve
4+
#
5+
# Import the Operating System so we can access the files for the VTS library
6+
from pythonnet import load
7+
load('coreclr')
8+
import clr
9+
import os
10+
file = '../libraries/Vts.dll'
11+
clr.AddReference(os.path.abspath(file))
12+
import numpy as np
13+
import plotly.graph_objects as go
14+
import plotly.express as px
15+
from Vts import *
16+
from Vts.Common import *
17+
from Vts.Extensions import *
18+
from Vts.Modeling.Optimizers import *
19+
from Vts.Modeling.ForwardSolvers import *
20+
from Vts.SpectralMapping import *
21+
from Vts.Factories import *
22+
from Vts.MonteCarlo import *
23+
from Vts.MonteCarlo.Sources import *
24+
from Vts.MonteCarlo.Tissues import *
25+
from Vts.MonteCarlo.Detectors import *
26+
from Vts.MonteCarlo.Factories import *
27+
from Vts.MonteCarlo.PhotonData import *
28+
from Vts.MonteCarlo.PostProcessing import *
29+
from System import Array, Object
30+
# Setup wavelengths in visible and NIR spectral regimes
31+
wavelengths = Array.CreateInstance(float, 8)
32+
for i in range(0, len(wavelengths)):
33+
wavelengths[i] = 650.0 + 50 * i
34+
# Setup fxs
35+
fxs = [0.0, 0.2]
36+
# Define parameters to fit [Hb,HbO2,A,b]
37+
measuredData = [28.4, 22.4, 1.2, 1.42]
38+
# Construct a scatter
39+
scattererMeasuredData = PowerLawScatterer(measuredData[2], measuredData[3])
40+
# Setup the values for the measured data
41+
chromophoresMeasuredData = Array.CreateInstance(IChromophoreAbsorber, 2)
42+
chromophoresMeasuredData[0] = ChromophoreAbsorber(ChromophoreType.HbO2, measuredData[0])
43+
chromophoresMeasuredData[1] = ChromophoreAbsorber(ChromophoreType.Hb, measuredData[1])
44+
opsMeasured = Tissue(chromophoresMeasuredData, scattererMeasuredData, "", n=1.4).GetOpticalProperties(wavelengths)
45+
# Create measurements using Nurbs-based white Monte Carlo forward solver
46+
measurementForwardSolver = NurbsForwardSolver()
47+
measuredROfFx= measurementForwardSolver.ROfFx(opsMeasured, fxs)
48+
# Create a forward solver as a model function for inversion
49+
forwardSolverForInversion = PointSourceSDAForwardSolver()
50+
51+
# Declare local forward reflectance function that computes reflectance
52+
# from chromophores and scatterer values
53+
# valuesSought = [Hb, HbO2, A, b]
54+
def CalculateReflectanceVsWavelengthFromChromophoreConcentration(
55+
valuesSought, wavelengths, fxs, forwardSolver):
56+
# Create a forward solve model function to solve inverse
57+
forwardSolverForInversion = PointSourceSDAForwardSolver()
58+
# Create an array of chromophore absorbers based on values
59+
chromophoresLocal = Array.CreateInstance(IChromophoreAbsorber, 2)
60+
chromophoresLocal[0] = ChromophoreAbsorber(ChromophoreType.HbO2, valuesSought[0])
61+
chromophoresLocal[1] = ChromophoreAbsorber(ChromophoreType.Hb, valuesSought[1])
62+
# Create a scatterr based on values
63+
scattererLocal = PowerLawScatterer(valuesSought[2], valuesSought[3])
64+
# Compose local tissue to obtain optical properties
65+
opsLocal = Tissue(chromophoresLocal, scattererLocal, "", n=1.4).GetOpticalProperties(wavelengths)
66+
print('opsLocal[0]=',opsLocal[0])
67+
# Compute reflectance for local absorbers
68+
modelDataLocal = Array.CreateInstance(float, len(wavelengths) * len(fxs))
69+
modelDataLocal = forwardSolver.ROfFx(opsLocal, fxs)
70+
print('modelDataLocal dims=',len(modelDataLocal))
71+
modelDataForReturn= Array.CreateInstance(float, len(wavelengths) * len(fxs))
72+
for i in range(0, len(wavelengths) * len(fxs)):
73+
modelDataForReturn[i] = modelDataLocal[i]
74+
return modelDataForReturn
75+
76+
# func for residual
77+
def residual(valuesSought, wavelengths, fxs, measuredROfFx, forwardSolver):
78+
predictedROfFx= CalculateReflectanceVsWavelengthFromChromophoreConcentration(
79+
valuesSought, wavelengths, fxs, forwardSolver)
80+
difference = Array.CreateInstance(float,len(wavelengths) * len(fxs))
81+
for i in range(0,len(wavelengths) * len(fxs)):
82+
difference[i] = predictedROfFx[i] - measuredROfFx[i]
83+
return difference
84+
85+
# Run the inversion: set up initial guess
86+
initialGuess = [18.0, 30.0, 0.8, 1.60]
87+
chromophoresInitialGuess = Array.CreateInstance(IChromophoreAbsorber, 2)
88+
chromophoresInitialGuess[0] = ChromophoreAbsorber(ChromophoreType.HbO2, initialGuess[0])
89+
chromophoresInitialGuess[1] = ChromophoreAbsorber(ChromophoreType.Hb, initialGuess[1])
90+
scattererInitialGuess = PowerLawScatterer(initialGuess[2], initialGuess[3])
91+
# Compose tissue for initial guess data to obtain OPs
92+
opsInitialGuess = Tissue(chromophoresInitialGuess, scattererInitialGuess, "", n=1.4).GetOpticalProperties(wavelengths)
93+
initialGuessROfFx = forwardSolverForInversion.ROfFx(opsInitialGuess, fxs)
94+
# Run the levenberg-marquardt inversion
95+
from scipy.optimize import least_squares
96+
fit = least_squares(
97+
residual,
98+
initialGuess,
99+
ftol=1e-9, xtol=1e-9, max_nfev=10000, # max_nfev needs to be integer
100+
args=(wavelengths, fxs, measuredROfFx, forwardSolverForInversion),
101+
method='lm')
102+
# Calculate final reflectance from model at fit values
103+
chromophoresFit = Array.CreateInstance(IChromophoreAbsorber, 2)
104+
chromophoresFit[0] = ChromophoreAbsorber(ChromophoreType.HbO2, fit.x[0])
105+
chromophoresFit[1] = ChromophoreAbsorber(ChromophoreType.Hb, fit.x[1])
106+
scattererFit = PowerLawScatterer(fit.x[2], fit.x[3])
107+
# Compose tissue for fit data to obtain OPs
108+
opsFit= Tissue(chromophoresFit, scattererFit, "", n=1.4).GetOpticalProperties(wavelengths)
109+
fitROfFx = forwardSolverForInversion.ROfFx(opsFit, fxs)
110+
# plot the results using Plotly
111+
xLabel = "wavelength [nm]"
112+
yLabel = "R(wavelength)"
113+
wvs = [w for w in wavelengths]
114+
# plot measured data
115+
meas = [m for m in measuredROfFx]
116+
chart = go.Figure()
117+
chart.add_trace(go.Scatter(x=wvs, y=meas, mode='markers', name='measured data'))
118+
# plot initial guess data
119+
ig = [i for i in initialGuessROfFx]
120+
chart.add_trace(go.Scatter(x=wvs, y=ig, mode='markers', name='initial guess'))
121+
# plot fit
122+
conv = [f for f in fitROfFx]
123+
chart.add_trace(go.Scatter(x=wvs, y=conv, mode='lines', name='converged'))
124+
chart.update_layout( title="ROfFx (inverse solution for chromophore concentrations, multiple wavelengths, multiple fx)", xaxis_title=xLabel, yaxis_title=yLabel)
125+
chart.show(renderer="browser")
126+
# output results
127+
print("Meas = [%5.3f %5.3f %5.3f %5.3f]" % (measuredData[0], measuredData[1], measuredData[2], measuredData[3]))
128+
print("IG = [%5.3f %5.3f %5.3f %5.3f] Chi2=%5.3e" % (
129+
initialGuess[0], initialGuess[1], initialGuess[2], initialGuess[3],
130+
np.dot(measuredROfFx,initialGuessROfFx)))
131+
print("Conv = [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (fit.x[0], fit.x[1], fit.x[2], fit.x[3],
132+
np.dot(measuredROfFx,fitROfFx)))
133+
print("error = [%5.3f %5.3f %5.3f %5.3f]" % (
134+
abs((measuredData[0]-fit.x[0])/measuredData[0]),
135+
abs((measuredData[1]-fit.x[1])/measuredData[1]),
136+
abs((measuredData[2]-fit.x[2])/measuredData[2]),
137+
abs((measuredData[3]-fit.x[3])/measuredData[3])))

0 commit comments

Comments
 (0)