Found error in Chi2 calculation and updated in all files. Modified iteration values to be the chromophores not mua.

Author: hayakawa16

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

@@ -67,6 +67,8 @@ def CalculateReflectanceVsWavelengthFromChromophoreConcAndScatterer(
    scattererLocal = PowerLawScatterer(valuesSought[2], valuesSought[3])
    # Compose local tissue to obtain optical properties
    opsLocal = Tissue(chromophoresLocal, scattererLocal, "", n=1.4).GetOpticalProperties(wavelengths)
+   print("iter:[Hb HbO2 A b]=[%5.3f %5.3f %5.3f %5.3f]" % (
+         valuesSought[0], valuesSought[1], valuesSought[2],valuesSought[3]))
    # Compute reflectance for local absorbers
    modelDataLocal=np.concatenate(
         [np.array(forwardSolver.ROfFx(opsLocal, fxs[0]), dtype=float), 
@@ -136,14 +138,15 @@ def residual(valuesSought, wavelengths, fxs, measuredROfFx, forwardSolver):
 
 chart.show(renderer="browser")
 # output results
-print("Meas =    [%5.3f %5.3f %5.3f %5.3f]" % (measuredData[0], measuredData[1], measuredData[2], measuredData[3]))
+print("Meas =    [%5.3f %5.3f %5.3f %5.3f]" % (
+                 measuredData[0], measuredData[1], measuredData[2], measuredData[3]))
 print("IG   =    [%5.3f %5.3f %5.3f %5.3f] Chi2=%5.3e" % (
                  initialGuess[0], initialGuess[1], initialGuess[2], initialGuess[3],
-                 np.dot(measuredROfFx,initialGuessROfFx)))
+                 np.dot(measuredROfFx-initialGuessROfFx,measuredROfFx-initialGuessROfFx)))
 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],
-                 np.dot(measuredROfFx,fitROfFx)))
-print("error =   [%5.3f %5.3f %5.3f %5.3f]" % (
-                 abs((measuredData[0]-fit.x[0])/measuredData[0]),
-                 abs((measuredData[1]-fit.x[1])/measuredData[1]),
-                 abs((measuredData[2]-fit.x[2])/measuredData[2]),
-                 abs((measuredData[3]-fit.x[3])/measuredData[3])))
+                 np.dot(measuredROfFx-fitROfFx,measuredROfFx-fitROfFx)))
+print("error =   [%3.2f %3.2f %3.2f %3.2f]%%" % (
+                 100*abs((measuredData[0]-fit.x[0])/measuredData[0]),
+                 100*abs((measuredData[1]-fit.x[1])/measuredData[1]),
+                 100*abs((measuredData[2]-fit.x[2])/measuredData[2]),
+                 100*abs((measuredData[3]-fit.x[3])/measuredData[3])))

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

@@ -60,7 +60,8 @@ def CalculateReflectanceFuncVsWavelengthFromChromophoreConcentration(
     chromophoresLocal[2] = ChromophoreAbsorber(ChromophoreType.H2O, chromophoreConcentration[2])
     # Compose local tissue to obtain optical properties
     opsLocal = Tissue(chromophoresLocal, params[2], "", n=1.4).GetOpticalProperties(params[0])
-    print('opsLocal[0]=',opsLocal[0])
+    print("iter:[HbO2,Hb,H2O]=[%3.2f %3.2f %3.2f]" % (
+       chromophoreConcentration[0],chromophoreConcentration[1],chromophoreConcentration[2]))
     # Compute reflectance for local absorbers
     modelDataLocal = forwardSolverForInversion.ROfRho(opsLocal, params[1]) 
     modelDataLocalCSharp = Array.CreateInstance(float, len(wavelengths))
@@ -120,11 +121,15 @@ def CalculateReflectanceFuncVsWavelengthFromChromophoreConcentration(
 chart.show(renderer="browser")
 # output results
 print("Meas =    [%5.3f %5.3f %5.3f]" % (measuredData[0], measuredData[1], measuredData[2]))
-print("IG   =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (initialGuess[0], initialGuess[1], initialGuess[2],
-                np.dot(measuredROfRho,initialGuessROfRho)))
+print("IG   =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (
+                initialGuess[0], initialGuess[1], initialGuess[2],
+                np.dot(np.subtract(measuredROfRho,initialGuessROfRho),
+                       np.subtract(measuredROfRho,initialGuessROfRho))))
 print("Conv =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (fit[0], fit[1], fit[2],
-                np.dot(measuredROfRho,fitROfRho)))
-print("error =   [%5.3f %5.3f %5.3f]" % (abs((measuredData[0]-fit[0])/measuredData[0]),
-                abs((measuredData[1]-fit[1])/measuredData[1]),
-                abs((measuredData[2]-fit[2])/measuredData[2])))
+                np.dot(np.subtract(measuredROfRho,fitROfRho),
+                       np.subtract(measuredROfRho,fitROfRho))))
+print("error =   [%5.3f %5.3f %5.3f]%%" % (
+                100*abs((measuredData[0]-fit[0])/measuredData[0]),
+                100*abs((measuredData[1]-fit[1])/measuredData[1]),
+                100*abs((measuredData[2]-fit[2])/measuredData[2])))
 

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

@@ -60,7 +60,8 @@ def CalculateReflectanceVsWavelengthFromChromophoreConcentration(
    chromophoresLocal[2] = ChromophoreAbsorber(ChromophoreType.H2O, chromophoreConcentration[2])
    # Compose local tissue to obtain optical properties
    opsLocal = Tissue(chromophoresLocal, scatterer, "", n=1.4).GetOpticalProperties(wavelengths)
-   print('opsLocal[0]=',opsLocal[0])
+   print("iter:[HbO2,Hb,H2O]=[%3.2f %3.2f %3.2f]" % (
+       chromophoreConcentration[0],chromophoreConcentration[1],chromophoreConcentration[2]))
    # Compute reflectance for local absorbers
    modelDataLocal = Array.CreateInstance(float, len(wavelengths))
    modelDataLocal = forwardSolver.ROfRho(opsLocal, rho) 
@@ -121,9 +122,12 @@ def residual(chromophoreConcentration, wavelengths, rho, scatterer, measuredROfR
 # output results
 print("Meas =    [%5.3f %5.3f %5.3f]" % (measuredData[0], measuredData[1], measuredData[2]))
 print("IG   =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (initialGuess[0], initialGuess[1], initialGuess[2],
-                np.dot(measuredROfRho,initialGuessROfRho)))
+                np.dot(np.subtract(measuredROfRho,initialGuessROfRho),
+                       np.subtract(measuredROfRho,initialGuessROfRho))))
 print("Conv =    [%5.3f %5.3f %5.3f] Chi2=%5.3e" % (fit.x[0], fit.x[1], fit.x[2],
-                np.dot(measuredROfRho,fitROfRho)))
-print("error =   [%5.3f %5.3f %5.3f]" % (abs((measuredData[0]-fit.x[0])/measuredData[0]),
-                abs((measuredData[1]-fit.x[1])/measuredData[1]),
-                abs((measuredData[2]-fit.x[2])/measuredData[2])))
+                np.dot(np.subtract(measuredROfRho,fitROfRho),
+                       np.subtract(measuredROfRho,fitROfRho))))
+print("error =   [%5.3f %5.3f %5.3f]%%" % (
+                100*abs((measuredData[0]-fit.x[0])/measuredData[0]),
+                100*abs((measuredData[1]-fit.x[1])/measuredData[1]),
+                100*abs((measuredData[2]-fit.x[2])/measuredData[2])))