From 6d944a7307fdb6d1e85b5db4ab57e7862761c167 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Tue, 1 Jul 2025 15:51:12 -0700 Subject: [PATCH 01/13] I took a whack at filling out the classes that need access to the newly added lamba0 parameter only needed by PowerLawScatterer (other IScatterers don't need). If this looks good, I will add unit tests. --- .../Spectroscopy/PowerLawScatterer.cs | 38 +++++++++++- src/Vts/Modeling/Spectroscopy/Tissue.cs | 58 ++++++++++++++++--- 2 files changed, 86 insertions(+), 10 deletions(-) diff --git a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs index e9c56dfd1..a49d1bbab 100644 --- a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs +++ b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs @@ -4,9 +4,10 @@ namespace Vts.SpectralMapping { /// /// Returns scattering values based on Steve Jacques' Skin Optics Summary: - /// https://omlc.ogi.edu/news/jan98/skinoptics.html + /// https://omlc.org/news/jan98/skinoptics.html /// This returned reduced scattering follows the approximate formula: - /// mus' = A1*lamda(-b1) + A2*lambda(-b2) + /// mus' = A1*(lamda/lambda0)(-b1) + A2*(lambda/lambda0)(-b2) + /// with default value of lambda0=1000nm /// public class PowerLawScatterer : BindableObject, IScatterer { @@ -14,6 +15,7 @@ public class PowerLawScatterer : BindableObject, IScatterer private double _b; private double _c; private double _d; + private double _lambda0; /// /// Constructs a power law scatterer; i.e. mus' = a*lamda^-b + c*lambda^-d @@ -28,6 +30,7 @@ public PowerLawScatterer(double a, double b, double c, double d) B = b; C = c; D = d; + Lambda0 = 1000; // use OMLC default value } /// @@ -40,6 +43,23 @@ public PowerLawScatterer(double a, double b) { } + /// + /// Constructs a power law scatterer; i.e. mus' = a*(lamda/lambda0)^-b + c*(lambda/lambda0)^-d + /// + /// The first prefactor + /// The first exponent + /// The second prefactor + /// The second exponent + /// Wavelength normalization factor + public PowerLawScatterer(double a, double b, double c, double d, double lambda0) + { + A = a; + B = b; + C = c; + D = d; + Lambda0 = lambda0; + } + /// /// Creates a power law scatterer using the specified tissue type /// @@ -70,6 +90,7 @@ public PowerLawScatterer() /// Tissue type public void SetTissueType(TissueType tissueType) { + Lambda0 = 1000; switch (tissueType) { case TissueType.Skin: @@ -176,6 +197,19 @@ public double D } } + /// + /// The wavelength normalization factor + /// + public double Lambda0 + { + get => _lambda0; + set + { + _lambda0 = value; + OnPropertyChanged("Lambda0"); + } + } + /// /// Returns mus' based on Steve Jacques' Skin Optics Summary: /// https://omlc.ogi.edu/news/jan98/skinoptics.html diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index 5335bd58d..985ea5ad0 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -93,12 +93,7 @@ public override string ToString() /// The absorption coefficient Mua public double GetMua(double wavelength) { - double mua = 0.0; - for (int i = 0; i < Absorbers.Count; i++) - { - mua += Absorbers[i].GetMua(wavelength); - } - return mua; + return Absorbers.Sum(t => t.GetMua(wavelength)); } /// @@ -111,6 +106,21 @@ public double GetMusp(double wavelength) return Scatterer != null ? Scatterer.GetMusp(wavelength) : 0; } + /// + /// Returns the reduced scattering coefficient for a given wavelength + /// + /// Wavelength + /// Wavelength normalization factor + /// The reduced scattering coefficient Mus' + public double GetMusp(double wavelength, double lambda0) + { + if (Scatterer != null) + { + return Scatterer is PowerLawScatterer scatterer ? scatterer.GetMusp(wavelength, lambda0) : Scatterer.GetMusp(wavelength); + } + return 0; + } + /// /// Returns the anisotropy coefficient for a given wavelength /// @@ -139,7 +149,7 @@ public double GetMus(double wavelength) public OpticalProperties GetOpticalProperties(double wavelength) { var mua = GetMua(wavelength); - var musp = GetMusp(wavelength); + var musp = GetMusp(wavelength, 1000); var g = GetG(wavelength); var n = N; return new OpticalProperties(mua, musp, g, n); @@ -148,18 +158,50 @@ public OpticalProperties GetOpticalProperties(double wavelength) /// /// Returns the optical properties for a given wavelength /// + /// Wavelength + /// Wavelength normalization factor + /// The optical properties + public OpticalProperties GetOpticalProperties(double wavelength, double lambda0) + { + var mua = GetMua(wavelength); + var musp = GetMusp(wavelength, lambda0); + var g = GetG(wavelength); + var n = N; + return new OpticalProperties(mua, musp, g, n); + } + + /// + /// Returns the optical properties for an array of wavelengths + /// /// Wavelength /// The optical properties public OpticalProperties[] GetOpticalProperties(double[] wavelengths) { var opArray = new OpticalProperties[wavelengths.Length]; - for (int i = 0; i < wavelengths.Length; i++) + for (var i = 0; i < wavelengths.Length; i++) { opArray[i] = GetOpticalProperties(wavelengths[i]); } return opArray; } + + /// + /// Returns the optical properties for an array of wavelengths + /// + /// Wavelength + /// Wavelength normalization factor + /// The optical properties + public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0) + { + var opArray = new OpticalProperties[wavelengths.Length]; + for (var i = 0; i < wavelengths.Length; i++) + { + opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); + } + return opArray; + } } + /// /// tissue provider class /// From 77aca78ad279a3dab0d0dc7a1ad0d82ee6551a25 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Tue, 1 Jul 2025 16:15:06 -0700 Subject: [PATCH 02/13] I had edited GetOpticalProperties without lambda0 parameter and should have. All existing unit tests pass. --- src/Vts/Modeling/Spectroscopy/Tissue.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index 985ea5ad0..b51c39a7e 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -149,7 +149,7 @@ public double GetMus(double wavelength) public OpticalProperties GetOpticalProperties(double wavelength) { var mua = GetMua(wavelength); - var musp = GetMusp(wavelength, 1000); + var musp = GetMusp(wavelength); var g = GetG(wavelength); var n = N; return new OpticalProperties(mua, musp, g, n); From 1a23e0346a12b52929ddfa506e130d20c19710e5 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Wed, 2 Jul 2025 18:03:22 -0700 Subject: [PATCH 03/13] Added unit tests for newly added code. --- .../Spectroscopy/PowerLawScattererTests.cs | 25 ++++++- .../Modeling/Spectroscopy/TissueTests.cs | 75 ++++++++++++++++++- 2 files changed, 98 insertions(+), 2 deletions(-) diff --git a/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs b/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs index fca648a58..9c7052520 100644 --- a/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs +++ b/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs @@ -7,6 +7,9 @@ namespace Vts.Test.Modeling.Spectroscopy [TestFixture] public class PowerLawScattererTests { + /// + /// Test constructor that specifies power law coefficients A, B, C, D + /// [Test] public void Test_power_law_scatterer_constructor() { @@ -18,6 +21,24 @@ public void Test_power_law_scatterer_constructor() Assert.That(scatterer.D, Is.EqualTo(0.0)); } + /// + /// Test constructor that specifies power law coefficients A, B, C, D, + /// and lambda0 the wavelength normalization factor + /// + [Test] + public void Test_power_law_scatterer_constructor_with_lambda0_specification() + { + var scatterer = new PowerLawScatterer(1, 0.1, 0.0, 0.0, 1000); + Assert.That(scatterer, Is.InstanceOf()); + Assert.That(scatterer.A, Is.EqualTo(1)); + Assert.That(scatterer.B, Is.EqualTo(0.1)); + Assert.That(scatterer.C, Is.EqualTo(0.0)); + Assert.That(scatterer.D, Is.EqualTo(0.0)); + } + + /// + /// Test that wrong specification of Tissue Type throws exception + /// [Test] public void Test_set_tissue_type_undefined() { @@ -25,6 +46,9 @@ public void Test_set_tissue_type_undefined() Assert.Throws(() => scatterer.SetTissueType((TissueType) 100)); } + /// + /// Test ability to call GetMusp method with and without lambda0 + /// [Test] public void Verify_user_ability_to_specify_lambda0() { @@ -38,7 +62,6 @@ public void Verify_user_ability_to_specify_lambda0() // set up call to GetMusp with another lambda0 musp = scatterer.GetMusp(1000, 800); Assert.That(musp, Is.Not.EqualTo(4.0)); - } } } diff --git a/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs b/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs index 5fd112a21..28a8e4b65 100644 --- a/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs +++ b/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs @@ -1,5 +1,4 @@ using NUnit.Framework; -using Vts.MonteCarlo; using Vts.SpectralMapping; namespace Vts.Test.Modeling.Spectroscopy @@ -9,6 +8,9 @@ public class TissueTests { private Tissue _tissue; + /// + /// Instantiate Tissue with predefined chromophores and scatterer + /// [OneTimeSetUp] public void One_time_setup() { @@ -28,6 +30,9 @@ public void One_time_setup() n); } + /// + /// Test setup specification of chromophores and optical properties + /// [Test] public void Test_tissue_constructor() { @@ -38,6 +43,9 @@ public void Test_tissue_constructor() Assert.That(_tissue.Absorbers[3].Concentration, Is.EqualTo(0.87)); } + /// + /// Test instantiation of Tissue with TissueTypes + /// [Test] public void Test_tissue_constructor_tissue_type() { @@ -52,6 +60,9 @@ public void Test_tissue_constructor_tissue_type() Assert.That(tissue.ScattererType, Is.EqualTo(ScatteringType.PowerLaw)); } + /// + /// Test GetMua method + /// [Test] public void Test_get_mua() { @@ -59,6 +70,9 @@ public void Test_get_mua() Assert.That(mua, Is.EqualTo(0.067854).Within(0.000001)); } + /// + /// Test GetMusp method + /// [Test] public void Test_get_musp() { @@ -66,6 +80,9 @@ public void Test_get_musp() Assert.That(musp, Is.EqualTo(0.839999).Within(0.000001)); } + /// + /// Test GetMus method + /// [Test] public void Test_get_mus() { @@ -73,6 +90,9 @@ public void Test_get_mus() Assert.That(mus, Is.EqualTo(4.2).Within(0.000001)); } + /// + /// Test GetG method + /// [Test] public void Test_get_g() { @@ -80,12 +100,18 @@ public void Test_get_g() Assert.That(g, Is.EqualTo(0.8)); } + /// + /// Test ToString method + /// [Test] public void Test_to_string() { Assert.That(_tissue.ToString(), Is.EqualTo("test_tissue")); } + /// + /// Test GetOpticalProperties method with single wavelength + /// [Test] public void Test_get_optical_properties() { @@ -97,6 +123,25 @@ public void Test_get_optical_properties() Assert.That(opticalProperties.G, Is.EqualTo(0.8)); } + /// + /// Test GetOpticalProperties method with single wavelength and lambda0 specification. + /// Lambda0 is the wavelength normalization specified in PowerLawScatterer + /// + [Test] + public void Test_get_optical_properties_with_lambda0_specification() + { + const double lambda0 = 1000; + var opticalProperties = _tissue.GetOpticalProperties(1000, lambda0); + Assert.That(opticalProperties.N, Is.EqualTo(1.4)); + Assert.That(opticalProperties.Mua, Is.EqualTo(0.067854).Within(0.000001)); + Assert.That(opticalProperties.Mus, Is.EqualTo(4.2).Within(0.000001)); + Assert.That(opticalProperties.Musp, Is.EqualTo(0.84).Within(0.000001)); + Assert.That(opticalProperties.G, Is.EqualTo(0.8)); + } + + /// + /// Test GetOpticalProperties with array of wavelengths + /// [Test] public void Test_get_optical_properties_wavelength_array() { @@ -119,5 +164,33 @@ public void Test_get_optical_properties_wavelength_array() Assert.That(opticalPropertyArray[2].Musp, Is.EqualTo(0.84).Within(0.000001)); Assert.That(opticalPropertyArray[2].G, Is.EqualTo(0.8)); } + + + /// + /// Test GetOpticalProperties with array of wavelengths and lambda0 specification + /// + [Test] + public void Test_get_optical_properties_wavelength_array_and_lambda0_specification() + { + const double lambda0 = 1000; + var wavelengths = new double[] { 600, 700, 1000 }; + var opticalPropertyArray = _tissue.GetOpticalProperties( + wavelengths, lambda0); + Assert.That(opticalPropertyArray[0].N, Is.EqualTo(1.4)); + Assert.That(opticalPropertyArray[0].Mua, Is.EqualTo(0.314619).Within(0.000001)); + Assert.That(opticalPropertyArray[0].Mus, Is.EqualTo(5.562449).Within(0.000001)); + Assert.That(opticalPropertyArray[0].Musp, Is.EqualTo(1.112489).Within(0.000001)); + Assert.That(opticalPropertyArray[0].G, Is.EqualTo(0.8)); + Assert.That(opticalPropertyArray[1].N, Is.EqualTo(1.4)); + Assert.That(opticalPropertyArray[1].Mua, Is.EqualTo(0.036097).Within(0.000001)); + Assert.That(opticalPropertyArray[1].Mus, Is.EqualTo(5.110287).Within(0.000001)); + Assert.That(opticalPropertyArray[1].Musp, Is.EqualTo(1.022057).Within(0.000001)); + Assert.That(opticalPropertyArray[1].G, Is.EqualTo(0.8)); + Assert.That(opticalPropertyArray[2].N, Is.EqualTo(1.4)); + Assert.That(opticalPropertyArray[2].Mua, Is.EqualTo(0.067854).Within(0.000001)); + Assert.That(opticalPropertyArray[2].Mus, Is.EqualTo(4.2).Within(0.000001)); + Assert.That(opticalPropertyArray[2].Musp, Is.EqualTo(0.84).Within(0.000001)); + Assert.That(opticalPropertyArray[2].G, Is.EqualTo(0.8)); + } } } From d80aa193368c22a88bdf0edeaf34318d70fddc47 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Tue, 8 Jul 2025 15:53:10 -0700 Subject: [PATCH 04/13] In PowerLawScatterer, put constructor with most parameters at top and modified subsequent constructors to call top constructor with fixed values. Same for GetMusp method. In Tissue cleaned up code. --- .../Spectroscopy/PowerLawScatterer.cs | 30 ++++++++----------- src/Vts/Modeling/Spectroscopy/Tissue.cs | 13 ++++---- 2 files changed, 20 insertions(+), 23 deletions(-) diff --git a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs index a49d1bbab..58c621b97 100644 --- a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs +++ b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs @@ -18,46 +18,42 @@ public class PowerLawScatterer : BindableObject, IScatterer private double _lambda0; /// - /// Constructs a power law scatterer; i.e. mus' = a*lamda^-b + c*lambda^-d + /// Constructs a power law scatterer; i.e. mus' = a*(lamda/lambda0)^-b + c*(lambda/lambda0)^-d /// /// The first prefactor /// The first exponent /// The second prefactor /// The second exponent - public PowerLawScatterer(double a, double b, double c, double d) + /// Wavelength normalization factor + public PowerLawScatterer(double a, double b, double c, double d, double lambda0) { A = a; B = b; C = c; D = d; - Lambda0 = 1000; // use OMLC default value + Lambda0 = lambda0; } /// - /// Creates a power law scatterer; i.e. mus' = a*lambda^-b + /// Constructs a power law scatterer; i.e. mus' = a*lamda^-b + c*lambda^-d /// /// The first prefactor /// The first exponent - public PowerLawScatterer(double a, double b) - : this(a,b,0.0,0.0) + /// The second prefactor + /// The second exponent + public PowerLawScatterer(double a, double b, double c, double d) + : this(a, b, c, d, 1000.0) { } /// - /// Constructs a power law scatterer; i.e. mus' = a*(lamda/lambda0)^-b + c*(lambda/lambda0)^-d + /// Creates a power law scatterer; i.e. mus' = a*lambda^-b /// /// The first prefactor /// The first exponent - /// The second prefactor - /// The second exponent - /// Wavelength normalization factor - public PowerLawScatterer(double a, double b, double c, double d, double lambda0) + public PowerLawScatterer(double a, double b) + : this(a,b,0.0,0.0, 1000.0) { - A = a; - B = b; - C = c; - D = d; - Lambda0 = lambda0; } /// @@ -218,7 +214,7 @@ public double Lambda0 /// The reduced scattering coefficient Mus' public double GetMusp(double wavelength) { - return A * Math.Pow(wavelength/1000, - B) + C * Math.Pow(wavelength/1000, - D); + return GetMusp(wavelength, 1000.0); } /// diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index b51c39a7e..fe28a555c 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -103,7 +103,7 @@ public double GetMua(double wavelength) /// The reduced scattering coefficient Mus' public double GetMusp(double wavelength) { - return Scatterer != null ? Scatterer.GetMusp(wavelength) : 0; + return Scatterer?.GetMusp(wavelength) ?? 0; } /// @@ -114,11 +114,12 @@ public double GetMusp(double wavelength) /// The reduced scattering coefficient Mus' public double GetMusp(double wavelength, double lambda0) { - if (Scatterer != null) + if (Scatterer == null) return 0; + if (double.IsNaN(lambda0)) { - return Scatterer is PowerLawScatterer scatterer ? scatterer.GetMusp(wavelength, lambda0) : Scatterer.GetMusp(wavelength); + return GetMusp(wavelength); } - return 0; + return Scatterer is PowerLawScatterer scatterer ? scatterer.GetMusp(wavelength, lambda0) : Scatterer.GetMusp(wavelength); } /// @@ -128,7 +129,7 @@ public double GetMusp(double wavelength, double lambda0) /// The anisotropy coefficient g public double GetG(double wavelength) { - return Scatterer != null ? Scatterer.GetG(wavelength) : 0; + return Scatterer?.GetG(wavelength) ?? 0; } /// @@ -138,7 +139,7 @@ public double GetG(double wavelength) /// The scattering coefficient Mus public double GetMus(double wavelength) { - return Scatterer != null ? Scatterer.GetMus(wavelength) : 0; + return Scatterer?.GetMus(wavelength) ?? 0; } /// From 4d52c36ed88299dff5d58cbb7f4bfbccd55cda41 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Tue, 8 Jul 2025 16:16:38 -0700 Subject: [PATCH 05/13] Corrected new Issue fround by sonarcloud. --- src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs index 58c621b97..b1889121b 100644 --- a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs +++ b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs @@ -202,7 +202,7 @@ public double Lambda0 set { _lambda0 = value; - OnPropertyChanged("Lambda0"); + OnPropertyChanged(nameof(Lambda0)); } } From 16fda69393be516a00153cda15370434d2902667 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Wed, 9 Jul 2025 11:16:53 -0700 Subject: [PATCH 06/13] Combined methods by making lambda0 be an optional parameter. --- src/Vts/Modeling/Spectroscopy/Tissue.cs | 51 +++++++++---------------- 1 file changed, 17 insertions(+), 34 deletions(-) diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index fe28a555c..59b3c684b 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -146,26 +146,12 @@ public double GetMus(double wavelength) /// Returns the optical properties for a given wavelength /// /// Wavelength + /// Optional wavelength normalization factor /// The optical properties - public OpticalProperties GetOpticalProperties(double wavelength) + public OpticalProperties GetOpticalProperties(double wavelength, double lambda0 = double.NaN) { + var musp = double.IsNaN(lambda0) ? GetMusp(wavelength) : GetMusp(wavelength, lambda0); var mua = GetMua(wavelength); - var musp = GetMusp(wavelength); - var g = GetG(wavelength); - var n = N; - return new OpticalProperties(mua, musp, g, n); - } - - /// - /// Returns the optical properties for a given wavelength - /// - /// Wavelength - /// Wavelength normalization factor - /// The optical properties - public OpticalProperties GetOpticalProperties(double wavelength, double lambda0) - { - var mua = GetMua(wavelength); - var musp = GetMusp(wavelength, lambda0); var g = GetG(wavelength); var n = N; return new OpticalProperties(mua, musp, g, n); @@ -175,29 +161,26 @@ public OpticalProperties GetOpticalProperties(double wavelength, double lambda0) /// Returns the optical properties for an array of wavelengths /// /// Wavelength + /// Optional wavelength normalization factor /// The optical properties - public OpticalProperties[] GetOpticalProperties(double[] wavelengths) + public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0 = double.NaN) { var opArray = new OpticalProperties[wavelengths.Length]; - for (var i = 0; i < wavelengths.Length; i++) + if (double.IsNaN(lambda0)) { - opArray[i] = GetOpticalProperties(wavelengths[i]); - } - return opArray; - } + for (var i = 0; i < wavelengths.Length; i++) + { + opArray[i] = GetOpticalProperties(wavelengths[i]); + } - /// - /// Returns the optical properties for an array of wavelengths - /// - /// Wavelength - /// Wavelength normalization factor - /// The optical properties - public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0) - { - var opArray = new OpticalProperties[wavelengths.Length]; - for (var i = 0; i < wavelengths.Length; i++) + } + else { - opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); + for (var i = 0; i < wavelengths.Length; i++) + { + opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); + } + } return opArray; } From d754ed22b7d6bc52d50d397d471aa1cef65f6977 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Mon, 14 Jul 2025 17:27:00 -0700 Subject: [PATCH 07/13] Replaced code that had an "if" statement outside of two "for" loops with code that has a single "for" and the "if" statement inside. This: https://stackoverflow.com/questions/38741599/using-an-if-statement-inside-loop-vs-outside-the-loop states that branch prediction will usually make sure that there is no difference at all in these caess. And the code is more readable. --- src/Vts/Modeling/Spectroscopy/Tissue.cs | 33 +++++++++++++++---------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index 59b3c684b..9c197bba6 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -166,22 +166,29 @@ public OpticalProperties GetOpticalProperties(double wavelength, double lambda0 public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0 = double.NaN) { var opArray = new OpticalProperties[wavelengths.Length]; - if (double.IsNaN(lambda0)) - { - for (var i = 0; i < wavelengths.Length; i++) - { - opArray[i] = GetOpticalProperties(wavelengths[i]); - } + //if (double.IsNaN(lambda0)) + //{ + // for (var i = 0; i < wavelengths.Length; i++) + // { + // opArray[i] = GetOpticalProperties(wavelengths[i]); + // } - } - else - { - for (var i = 0; i < wavelengths.Length; i++) - { - opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); - } + //} + //else + //{ + // for (var i = 0; i < wavelengths.Length; i++) + // { + // opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); + // } + //} + for (var i = 0; i < wavelengths.Length; i++) + { + opArray[i] = double.IsNaN(lambda0) + ? GetOpticalProperties(wavelengths[i]) + : GetOpticalProperties(wavelengths[i], lambda0); } + return opArray; } } From b5524e1bad815daf71bfb5549af234e46e9804b5 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Mon, 14 Jul 2025 17:48:01 -0700 Subject: [PATCH 08/13] Revert "Replaced code that had an "if" statement outside of two "for" loops with code that has a single "for" and the "if" statement inside. This:" This reverts commit d754ed22b7d6bc52d50d397d471aa1cef65f6977. --- src/Vts/Modeling/Spectroscopy/Tissue.cs | 33 ++++++++++--------------- 1 file changed, 13 insertions(+), 20 deletions(-) diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index 9c197bba6..59b3c684b 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -166,29 +166,22 @@ public OpticalProperties GetOpticalProperties(double wavelength, double lambda0 public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0 = double.NaN) { var opArray = new OpticalProperties[wavelengths.Length]; - //if (double.IsNaN(lambda0)) - //{ - // for (var i = 0; i < wavelengths.Length; i++) - // { - // opArray[i] = GetOpticalProperties(wavelengths[i]); - // } - - //} - //else - //{ - // for (var i = 0; i < wavelengths.Length; i++) - // { - // opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); - // } - - //} - for (var i = 0; i < wavelengths.Length; i++) + if (double.IsNaN(lambda0)) { - opArray[i] = double.IsNaN(lambda0) - ? GetOpticalProperties(wavelengths[i]) - : GetOpticalProperties(wavelengths[i], lambda0); + for (var i = 0; i < wavelengths.Length; i++) + { + opArray[i] = GetOpticalProperties(wavelengths[i]); + } + } + else + { + for (var i = 0; i < wavelengths.Length; i++) + { + opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); + } + } return opArray; } } From d049a83fcbe65659d48ad1f1ee3d223c68673fab Mon Sep 17 00:00:00 2001 From: hayakawa Date: Mon, 14 Jul 2025 17:58:22 -0700 Subject: [PATCH 09/13] Replaced code that had an "if" statement outside of two "for" loops with code that has a single "for" and the "if" statement inside. This: https://stackoverflow.com/questions/38741599/using-an-if-statement-inside-loop-vs-outside-the-loop states that branch prediction will usually make sure that there is no difference at all in these caess. And the code is more readable. --- src/Vts/Modeling/Spectroscopy/Tissue.cs | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index 59b3c684b..7cebea712 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -166,21 +166,11 @@ public OpticalProperties GetOpticalProperties(double wavelength, double lambda0 public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0 = double.NaN) { var opArray = new OpticalProperties[wavelengths.Length]; - if (double.IsNaN(lambda0)) + for (var i = 0; i < wavelengths.Length; i++) { - for (var i = 0; i < wavelengths.Length; i++) - { - opArray[i] = GetOpticalProperties(wavelengths[i]); - } - - } - else - { - for (var i = 0; i < wavelengths.Length; i++) - { - opArray[i] = GetOpticalProperties(wavelengths[i], lambda0); - } - + opArray[i] = double.IsNaN(lambda0) + ? GetOpticalProperties(wavelengths[i]) + : GetOpticalProperties(wavelengths[i], lambda0); } return opArray; } From 1bb93004df8ef8788de7458858f57d31bf913c02 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 7 Aug 2025 16:35:52 -0700 Subject: [PATCH 10/13] Removed overloads and injected lambda0 via the instantiation of PowerLawScatterer. I'm not sure why I didn't see this solution earlier. Updated unit tests that test using default lambda0 value (1000nm) and user specified value. --- .../Spectroscopy/PowerLawScattererTests.cs | 8 +-- .../Modeling/Spectroscopy/TissueTests.cs | 63 ++++++++++++------- .../Spectroscopy/PowerLawScatterer.cs | 14 +---- src/Vts/Modeling/Spectroscopy/Tissue.cs | 28 ++------- 4 files changed, 51 insertions(+), 62 deletions(-) diff --git a/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs b/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs index 9c7052520..982d13237 100644 --- a/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs +++ b/src/Vts.Test/Modeling/Spectroscopy/PowerLawScattererTests.cs @@ -53,14 +53,12 @@ public void Test_set_tissue_type_undefined() public void Verify_user_ability_to_specify_lambda0() { // set up call to GetMusp without lambda0 specified - var scatterer = new PowerLawScatterer(1.0, 2.0, 3.0, 4.0); + var scatterer = new PowerLawScatterer(1.0, 2.0, 3.0, 4.0, 1000.0); var musp = scatterer.GetMusp(1000); Assert.That(musp, Is.EqualTo(4.0)); - // set up call to GetMusp with lambda0 = 1000 specified - musp = scatterer.GetMusp(1000, 1000); - Assert.That(musp, Is.EqualTo(4.0)); // set up call to GetMusp with another lambda0 - musp = scatterer.GetMusp(1000, 800); + scatterer = new PowerLawScatterer(1.0, 2.0, 3.0, 4.0, 800.0); + musp = scatterer.GetMusp(1000); Assert.That(musp, Is.Not.EqualTo(4.0)); } } diff --git a/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs b/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs index 28a8e4b65..d70d69aca 100644 --- a/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs +++ b/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs @@ -14,7 +14,7 @@ public class TissueTests [OneTimeSetUp] public void One_time_setup() { - // used values for tissue=liver + // used values for tissue=liver, power law wavelength normalization set to default = 1000nm var scatterer = new PowerLawScatterer(0.84, 0.55); var hbAbsorber = new ChromophoreAbsorber(ChromophoreType.Hb, 66); var hbo2Absorber = new ChromophoreAbsorber(ChromophoreType.HbO2, 124); @@ -124,18 +124,29 @@ public void Test_get_optical_properties() } /// - /// Test GetOpticalProperties method with single wavelength and lambda0 specification. - /// Lambda0 is the wavelength normalization specified in PowerLawScatterer + /// Test Tissue with PowerLawScatterer method with single wavelength and PowerLawScatterer + /// wavelength normalization lambda0 specified to non-default value /// [Test] - public void Test_get_optical_properties_with_lambda0_specification() + public void Test_power_law_scatterer_with_lambda0_specification() { - const double lambda0 = 1000; - var opticalProperties = _tissue.GetOpticalProperties(1000, lambda0); + // used values for tissue=liver, power law wavelength normalization set to default = 1000nm + var scatterer = new PowerLawScatterer(0.84, 0.55, 0.0, 0.0, 750.0); + var hbAbsorber = new ChromophoreAbsorber(ChromophoreType.Hb, 66); + var hbo2Absorber = new ChromophoreAbsorber(ChromophoreType.HbO2, 124); + var fatAbsorber = new ChromophoreAbsorber(ChromophoreType.Fat, 0.02); + var waterAbsorber = new ChromophoreAbsorber(ChromophoreType.H2O, 0.87); + const double n = 1.4; + _tissue = new Tissue( + new IChromophoreAbsorber[] { hbAbsorber, hbo2Absorber, fatAbsorber, waterAbsorber }, + scatterer, + "test_tissue", + n); + var opticalProperties = _tissue.GetOpticalProperties(1000); Assert.That(opticalProperties.N, Is.EqualTo(1.4)); Assert.That(opticalProperties.Mua, Is.EqualTo(0.067854).Within(0.000001)); - Assert.That(opticalProperties.Mus, Is.EqualTo(4.2).Within(0.000001)); - Assert.That(opticalProperties.Musp, Is.EqualTo(0.84).Within(0.000001)); + Assert.That(opticalProperties.Mus, Is.EqualTo(3.585361).Within(0.000001)); + Assert.That(opticalProperties.Musp, Is.EqualTo(0.717073).Within(0.000001)); Assert.That(opticalProperties.G, Is.EqualTo(0.8)); } @@ -165,31 +176,41 @@ public void Test_get_optical_properties_wavelength_array() Assert.That(opticalPropertyArray[2].G, Is.EqualTo(0.8)); } - /// - /// Test GetOpticalProperties with array of wavelengths and lambda0 specification + /// Test Tissue with PowerLawScatterer method with multiple wavelengths and PowerLawScatterer + /// wavelength normalization lambda0 specified to non-default value /// [Test] - public void Test_get_optical_properties_wavelength_array_and_lambda0_specification() - { - const double lambda0 = 1000; + public void Test_power_law_scatterer_with_lambda0_specification_array_of_wavelengths() + { + // used values for tissue=liver, power law wavelength normalization = 750nm + var scatterer = new PowerLawScatterer(0.84, 0.55, 0.0, 0.0, 750.0); + var hbAbsorber = new ChromophoreAbsorber(ChromophoreType.Hb, 66); + var hbo2Absorber = new ChromophoreAbsorber(ChromophoreType.HbO2, 124); + var fatAbsorber = new ChromophoreAbsorber(ChromophoreType.Fat, 0.02); + var waterAbsorber = new ChromophoreAbsorber(ChromophoreType.H2O, 0.87); + const double n = 1.4; + _tissue = new Tissue( + new IChromophoreAbsorber[] { hbAbsorber, hbo2Absorber, fatAbsorber, waterAbsorber }, + scatterer, + "test_tissue", + n); var wavelengths = new double[] { 600, 700, 1000 }; - var opticalPropertyArray = _tissue.GetOpticalProperties( - wavelengths, lambda0); + var opticalPropertyArray = _tissue.GetOpticalProperties(wavelengths); Assert.That(opticalPropertyArray[0].N, Is.EqualTo(1.4)); Assert.That(opticalPropertyArray[0].Mua, Is.EqualTo(0.314619).Within(0.000001)); - Assert.That(opticalPropertyArray[0].Mus, Is.EqualTo(5.562449).Within(0.000001)); - Assert.That(opticalPropertyArray[0].Musp, Is.EqualTo(1.112489).Within(0.000001)); + Assert.That(opticalPropertyArray[0].Mus, Is.EqualTo(4.748427).Within(0.000001)); + Assert.That(opticalPropertyArray[0].Musp, Is.EqualTo(0.949685).Within(0.000001)); Assert.That(opticalPropertyArray[0].G, Is.EqualTo(0.8)); Assert.That(opticalPropertyArray[1].N, Is.EqualTo(1.4)); Assert.That(opticalPropertyArray[1].Mua, Is.EqualTo(0.036097).Within(0.000001)); - Assert.That(opticalPropertyArray[1].Mus, Is.EqualTo(5.110287).Within(0.000001)); - Assert.That(opticalPropertyArray[1].Musp, Is.EqualTo(1.022057).Within(0.000001)); + Assert.That(opticalPropertyArray[1].Mus, Is.EqualTo(4.362435).Within(0.000001)); + Assert.That(opticalPropertyArray[1].Musp, Is.EqualTo(0.872487).Within(0.000001)); Assert.That(opticalPropertyArray[1].G, Is.EqualTo(0.8)); Assert.That(opticalPropertyArray[2].N, Is.EqualTo(1.4)); Assert.That(opticalPropertyArray[2].Mua, Is.EqualTo(0.067854).Within(0.000001)); - Assert.That(opticalPropertyArray[2].Mus, Is.EqualTo(4.2).Within(0.000001)); - Assert.That(opticalPropertyArray[2].Musp, Is.EqualTo(0.84).Within(0.000001)); + Assert.That(opticalPropertyArray[2].Mus, Is.EqualTo(3.585361).Within(0.000001)); + Assert.That(opticalPropertyArray[2].Musp, Is.EqualTo(0.717072).Within(0.000001)); Assert.That(opticalPropertyArray[2].G, Is.EqualTo(0.8)); } } diff --git a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs index b1889121b..d2cdbd6ae 100644 --- a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs +++ b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs @@ -209,23 +209,13 @@ public double Lambda0 /// /// Returns mus' based on Steve Jacques' Skin Optics Summary: /// https://omlc.ogi.edu/news/jan98/skinoptics.html + /// and normalizes wavelength by _lambda0 (default=1000[nm]) /// /// Wavelength /// The reduced scattering coefficient Mus' public double GetMusp(double wavelength) { - return GetMusp(wavelength, 1000.0); - } - - /// - /// Returns mus' given wavelength and lambda0 specification - /// - /// Wavelength - /// Wavelength normalization factor - /// The reduced scattering coefficient Mus' - public double GetMusp(double wavelength, double lambda0) - { - return A * Math.Pow(wavelength / lambda0, -B) + C * Math.Pow(wavelength / lambda0, -D); + return A * Math.Pow(wavelength / _lambda0, -B) + C * Math.Pow(wavelength / _lambda0, -D); } /// diff --git a/src/Vts/Modeling/Spectroscopy/Tissue.cs b/src/Vts/Modeling/Spectroscopy/Tissue.cs index 7cebea712..04322c7f0 100644 --- a/src/Vts/Modeling/Spectroscopy/Tissue.cs +++ b/src/Vts/Modeling/Spectroscopy/Tissue.cs @@ -106,22 +106,6 @@ public double GetMusp(double wavelength) return Scatterer?.GetMusp(wavelength) ?? 0; } - /// - /// Returns the reduced scattering coefficient for a given wavelength - /// - /// Wavelength - /// Wavelength normalization factor - /// The reduced scattering coefficient Mus' - public double GetMusp(double wavelength, double lambda0) - { - if (Scatterer == null) return 0; - if (double.IsNaN(lambda0)) - { - return GetMusp(wavelength); - } - return Scatterer is PowerLawScatterer scatterer ? scatterer.GetMusp(wavelength, lambda0) : Scatterer.GetMusp(wavelength); - } - /// /// Returns the anisotropy coefficient for a given wavelength /// @@ -146,11 +130,10 @@ public double GetMus(double wavelength) /// Returns the optical properties for a given wavelength /// /// Wavelength - /// Optional wavelength normalization factor /// The optical properties - public OpticalProperties GetOpticalProperties(double wavelength, double lambda0 = double.NaN) + public OpticalProperties GetOpticalProperties(double wavelength) { - var musp = double.IsNaN(lambda0) ? GetMusp(wavelength) : GetMusp(wavelength, lambda0); + var musp = GetMusp(wavelength); var mua = GetMua(wavelength); var g = GetG(wavelength); var n = N; @@ -161,16 +144,13 @@ public OpticalProperties GetOpticalProperties(double wavelength, double lambda0 /// Returns the optical properties for an array of wavelengths /// /// Wavelength - /// Optional wavelength normalization factor /// The optical properties - public OpticalProperties[] GetOpticalProperties(double[] wavelengths, double lambda0 = double.NaN) + public OpticalProperties[] GetOpticalProperties(double[] wavelengths) { var opArray = new OpticalProperties[wavelengths.Length]; for (var i = 0; i < wavelengths.Length; i++) { - opArray[i] = double.IsNaN(lambda0) - ? GetOpticalProperties(wavelengths[i]) - : GetOpticalProperties(wavelengths[i], lambda0); + opArray[i] = GetOpticalProperties(wavelengths[i]); } return opArray; } From 9a3521b553912d163d687761eeaf8ab5d00bd1fe Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 7 Aug 2025 16:45:52 -0700 Subject: [PATCH 11/13] Added overload to PowerLawScatterer that allows user to specify A(lambda/lambda0)^-b, i.e. simple not extended power law definition, --- src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs index d2cdbd6ae..74ea85c70 100644 --- a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs +++ b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs @@ -46,6 +46,17 @@ public PowerLawScatterer(double a, double b, double c, double d) { } + /// + /// Creates a power law scatterer; i.e. mus' = a*(lambda/lambda0)^-b + /// + /// The first prefactor + /// The first exponent + /// Wavelength normalization factor + public PowerLawScatterer(double a, double b, double lambda0) + : this(a, b, 0.0, 0.0, lambda0) + { + } + /// /// Creates a power law scatterer; i.e. mus' = a*lambda^-b /// From 7158b8d211b17090bdd3da478b110ac6550c37ce Mon Sep 17 00:00:00 2001 From: hayakawa Date: Thu, 7 Aug 2025 16:54:47 -0700 Subject: [PATCH 12/13] Modified test to use new PowerLawScatterer constructor. --- src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs b/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs index d70d69aca..b91287ce0 100644 --- a/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs +++ b/src/Vts.Test/Modeling/Spectroscopy/TissueTests.cs @@ -184,7 +184,7 @@ public void Test_get_optical_properties_wavelength_array() public void Test_power_law_scatterer_with_lambda0_specification_array_of_wavelengths() { // used values for tissue=liver, power law wavelength normalization = 750nm - var scatterer = new PowerLawScatterer(0.84, 0.55, 0.0, 0.0, 750.0); + var scatterer = new PowerLawScatterer(0.84, 0.55, 750.0); var hbAbsorber = new ChromophoreAbsorber(ChromophoreType.Hb, 66); var hbo2Absorber = new ChromophoreAbsorber(ChromophoreType.HbO2, 124); var fatAbsorber = new ChromophoreAbsorber(ChromophoreType.Fat, 0.02); From cc1b454e750b281dc4010bb651c548a64dff7f19 Mon Sep 17 00:00:00 2001 From: hayakawa Date: Sun, 10 Aug 2025 16:39:12 -0700 Subject: [PATCH 13/13] Combined constructor overloads down to 2. There are two use cases I needed to cover: 1. mus' = A1*(lamda/lambda0)(-b1) + A2*(lambda/lambda0)(-b2) with an optional lambda0, i.e. public PowerLawScatterer(double a, double b, double c, double d, double lambda0 = 1000.0). 2. mus' = a*(lamda/lambda0)^-b with an optional lambda0, i.e. public PowerLawScatterer(double a, double b, double lambda0 = 1000.0). I couldn't combine without parameter ambiguity. --- .../Spectroscopy/PowerLawScatterer.cs | 28 ++----------------- 1 file changed, 3 insertions(+), 25 deletions(-) diff --git a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs index 74ea85c70..8171aad61 100644 --- a/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs +++ b/src/Vts/Modeling/Spectroscopy/PowerLawScatterer.cs @@ -25,7 +25,7 @@ public class PowerLawScatterer : BindableObject, IScatterer /// The second prefactor /// The second exponent /// Wavelength normalization factor - public PowerLawScatterer(double a, double b, double c, double d, double lambda0) + public PowerLawScatterer(double a, double b, double c, double d, double lambda0 = 1000.0) { A = a; B = b; @@ -35,38 +35,16 @@ public PowerLawScatterer(double a, double b, double c, double d, double lambda0) } /// - /// Constructs a power law scatterer; i.e. mus' = a*lamda^-b + c*lambda^-d - /// - /// The first prefactor - /// The first exponent - /// The second prefactor - /// The second exponent - public PowerLawScatterer(double a, double b, double c, double d) - : this(a, b, c, d, 1000.0) - { - } - - /// - /// Creates a power law scatterer; i.e. mus' = a*(lambda/lambda0)^-b + /// Constructs a power law scatterer; i.e. mus' = a*(lamda/lambda0)^-b /// /// The first prefactor /// The first exponent /// Wavelength normalization factor - public PowerLawScatterer(double a, double b, double lambda0) + public PowerLawScatterer(double a, double b, double lambda0 = 1000.0) : this(a, b, 0.0, 0.0, lambda0) { } - /// - /// Creates a power law scatterer; i.e. mus' = a*lambda^-b - /// - /// The first prefactor - /// The first exponent - public PowerLawScatterer(double a, double b) - : this(a,b,0.0,0.0, 1000.0) - { - } - /// /// Creates a power law scatterer using the specified tissue type ///