diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 000000000..90d995bff --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "dotnet.defaultSolution": "src/Vts.sln" +} \ No newline at end of file diff --git a/matlab/vts_wrapper/short_course_monte_carlo_lab.m b/matlab/vts_wrapper/short_course_monte_carlo_lab.m index 2f48cfc4a..d3d2ff5cf 100644 --- a/matlab/vts_wrapper/short_course_monte_carlo_lab.m +++ b/matlab/vts_wrapper/short_course_monte_carlo_lab.m @@ -7,7 +7,7 @@ startup(); % ======================================================================= % -%% Example 1a: run Monte Carlo simulations for fluence with increasing N photons +%% Example 01a: run Monte Carlo simulations for fluence with increasing N photons clear; Nphot=[10, 100, 1000, 10000]; % number of photons launched takes about 1 mins @@ -129,7 +129,8 @@ save(sprintf('dawSD%d.txt',j),'-ascii','SD'); end end -%% Example 1b: Compare standard deviation of two absorption methods + +%% Example 01b: Compare standard deviation of two absorption methods % this cell relies on above cell execution f=figure('Position',[scrsz(3)/5 scrsz(4)/5 scrsz(3)/2 scrsz(4)/2]); set(f,'Name','Analog relative error - DAW relative error'); @@ -156,7 +157,7 @@ end % ======================================================================= % -%% Example 2: run Monte Carlo simulations accounting for absorption with +%% Example 02: run Monte Carlo simulations accounting for absorption with % analog and continuous absorption weighting with 10,000 photons and compare % time and relative error diff --git a/matlab/vts_wrapper/vts_mc_demo.m b/matlab/vts_wrapper/vts_mc_demo.m index 60289d67f..fc0332b97 100644 --- a/matlab/vts_wrapper/vts_mc_demo.m +++ b/matlab/vts_wrapper/vts_mc_demo.m @@ -8,7 +8,7 @@ startup(); % ======================================================================= % -%% Example 1: run a simple Monte Carlo simulation with 1000 photons +%% Example 01: run a simple Monte Carlo simulation with 1000 photons % create a default set of inputs si = SimulationInput(); @@ -29,7 +29,7 @@ set(f,'Name','Reflectance vs Rho for N=1000'); % ======================================================================= % -%% Example 2: run Monte Carlo simulations for two absorption weighting types +%% Example 02: run Monte Carlo simulations for two absorption weighting types % with 1000 photons each and compare computation time % create a default set of inputs @@ -49,7 +49,7 @@ output2 = VtsMonteCarlo.RunSimulation(si); % ======================================================================= % -%% Example 3: run a Monte Carlo simulation with a fully-customized input +%% Example 03: run a Monte Carlo simulation with a fully-customized input % (values used here are the class defaults) % 1) define a source... @@ -128,7 +128,7 @@ output = VtsMonteCarlo.RunSimulation(input); % ======================================================================= % -%% Example 4: run a list of Monte Carlo simulations +%% Example 04: run a list of Monte Carlo simulations % create a list of two default SimulationInput with different numbers of % photons @@ -156,7 +156,7 @@ set(f,'Name','Simulation 2 results'); % ======================================================================= % -% Example 5: run a Monte Carlo simulation with post-processing enabled +% Example 05: run a Monte Carlo simulation with post-processing enabled % First run a simulation, then post-process the generated database and % compare on-the-fly results with post-processed results si = SimulationInput(); @@ -187,7 +187,7 @@ legend('on-the-fly','post-processed'); % ======================================================================= % -% Example 6: run a Monte Carlo simulation with pMC post-processing enabled +% Example 06: run a Monte Carlo simulation with pMC post-processing enabled % First run a simulation, then post-process the generated database CA % varying optical properties si = SimulationInput(); @@ -260,7 +260,7 @@ set(f,'Name','Perturbation MC using CAW'); % ======================================================================= % -% Example 7: run a Monte Carlo simulation with pMC post-processing enabled +% Example 07: run a Monte Carlo simulation with pMC post-processing enabled % Use generated database to solve inverse problem with measured data % generated using Nurbs % NOTE: convergence to measured data optical properties affected by: @@ -375,7 +375,7 @@ abs(measMus-recoveredOPs(2))/measOPs(2))); % ======================================================================= % -% Example 8: run a Monte Carlo simulation and verify results with +% Example 08: run a Monte Carlo simulation and verify results with % those in unit tests in Visual Studio % spell out all input to ensure same settings as in unit test si = SimulationInput(); @@ -422,8 +422,9 @@ (abs(d3.Mean(1,1)-95.492965855)<0.000000001)) disp('unit tests pass'); end + % ======================================================================= % -%% Example 9: run a Monte Carlo simulation for transmittance tallies with 1000 photons +%% Example 09: run a Monte Carlo simulation for transmittance tallies with 1000 photons % create a default set of inputs si = SimulationInput(); @@ -463,6 +464,7 @@ f = figure; plot(d2.Fx, abs(d2.Mean)); ylabel('T(fx) [unitless]'); xlabel('Fx (mm^-^1)'); title('Transmittance vs fx for N=1000'); set(f,'Name','Transmittance vs Fx for N=1000'); + % ======================================================================= % %% Example 10: run R(fx) detector results diff --git a/matlab/vts_wrapper/vts_solver_demo.m b/matlab/vts_wrapper/vts_solver_demo.m index 9c2d27907..ad7c037ad 100644 --- a/matlab/vts_wrapper/vts_solver_demo.m +++ b/matlab/vts_wrapper/vts_solver_demo.m @@ -8,7 +8,7 @@ startup(); -%% Example ROfRhoAndFt +%% Example 01: ROfRhoAndFt % Evaluate reflectance as a function of rho and temporal-frequency with one % set of optical properites. @@ -41,7 +41,7 @@ xlabel('f_t'); set(f,'Name','Reflectance as a function of Rho and Temporal-frequency'); -%% Example ROfFxAndFt +%% Example 02: ROfFxAndFt % Evaluate reflectance as a function of spatial- and temporal- frequencies % with one set of optical properites. @@ -68,7 +68,7 @@ xlabel('f_t'); set(f,'Name','Reflectance as a function of Spatial- and Temporal- frequencies'); -%% Example FluenceOfRhoAndZ +%% Example 03: FluenceOfRhoAndZ % Evaluate fluence as a function of rho and z using optical properties from % a list of chromophore absorbers with their concentrations and a power law % scatterer for a range of wavelengths. @@ -105,11 +105,11 @@ f = figure; imagesc(wv,zs,log(squeeze(test(:,1,:)))); ylabel('z [mm]'); xlabel('wavelength, \lambda [nm]'); -title('Fluence of \lambda and z and \rho=0.1 mm'); -set(f,'Name','Fluence as a function of Rho and z'); +title('Fluence of \lambda and z at \rho=0.1 mm'); +set(f,'Name','Fluence as a function of lambda and z'); -%% Example FluenceOfRhoAndZAndFt +%% Example 04: FluenceOfRhoAndZAndFt % Evaluate fluence as a function of rho and z using one set of optical % properties and a distributed gaussian source SDA solver type. @@ -148,7 +148,7 @@ ylabel('z [mm]') set(f,'Name','Modulation of fluence (AC/DC) of Rho and z and ft (ft=1GHz)'); -%% Example PHDOfRhoAndZ +%% Example 05: PHDOfRhoAndZ % Evaluate Photon Hitting Density in cylindrical coordinates % using one set of optical properties and a distributed gaussian source SDA % solver type. @@ -167,7 +167,7 @@ ylabel('z [mm]') set(f,'Name','PHD of Rho and z'); -%% Example FluenceOfRhoAndZTwoLayer +%% Example 06: FluenceOfRhoAndZTwoLayer % Evaluate fluence in cylindrical coordinates for a two % layer tissue with specified source-detector separation and top layer thickness % using two sets of optical properties. @@ -185,7 +185,7 @@ ylabel('z [mm]'); set(f,'Name','Fluence of Rho and z - Two Layer'); -%% Example PHDOfRhoAndZTwoLayer +%% Example 07: PHDOfRhoAndZTwoLayer % Evaluate Photon Hitting Density in cylindrical coordinates for a two % layer tissue with specified source-detector separation and top layer thickness % using two sets of optical properties. @@ -204,7 +204,7 @@ ylabel('z [mm]'); set(f,'Name','PHD of Rho and z - Two Layer'); -%% Example AbsorbedEnergyOfRhoAndZ +%% Example 08: AbsorbedEnergyOfRhoAndZ % Evaluate Absorbed Energy of rho and z using one set of optical properties % and a point source SDA solver type. @@ -225,7 +225,7 @@ ylabel('z [mm]'); set(f,'Name','Absorbed Energy of Rho and z'); -%% Example ROfRho +%% Example 09: ROfRho % Evaluate reflectance as a function of rho with three sets % of optical properites. @@ -242,7 +242,7 @@ ylabel('R(\rho)'); xlabel('\rho'); -%% Example ROfRhoAndT +%% Example 10: ROfRhoAndTime % Evaluate reflectance of rho and t at one s-d separation and % two sets of optical properites. @@ -259,7 +259,7 @@ ylabel('R(t)'); xlabel('Time, t [ns]'); -%% Example ROfFxAndT +%% Example 11: ROfFxAndTime % Evaluate reflectance as a function of spacial-frequency and t using one % set of optical properites. @@ -291,7 +291,7 @@ end legend(l2, 'FontSize', 12); -%% Example ROfFx (single set of optical properties) +%% Example 12: ROfFx (single set of optical properties) % Evaluate reflectance as a function of spacial-frequency with a single set % of optical properties. @@ -304,7 +304,7 @@ ylabel('R(f_x)'); xlabel('Spatial frequency, f_x [mm^-^1]'); -%% Example ROfFx (multiple sets of optical properties) +%% Example 13: ROfFx (multiple sets of optical properties) % Evaluate reflectance as a function of spacial-frequency % with multiple sets of optical properties, varying mua linearly. @@ -323,8 +323,7 @@ ylabel('R(f_x)'); xlabel('Spatial frequency, f_x [mm^-^1]'); - -%% Example ROfFx (multiple optical properties, varying mua as a function of wavelength, mus' as intralipid scatterer) +%% Example 14: ROfFx (multiple optical properties, varying mua as a function of wavelength, mus' as intralipid scatterer) % Evaluate reflectance as a function of spacial-frequency with multiple % sets of optical properties, varying mua as a function of wavelength. @@ -343,14 +342,14 @@ % % or scatterer.Type = 'Intralipid'; -scatterer.vol_frac = 0.5; +scatterer.vol_frac = 0.01; % % or % scatterer.Type = 'Mie'; % scatterer.radius = 0.5; % scatterer.n = 1.4; % scatterer.nMedium = 1.0; -% scatterer.vol_frac = 0.5; +% scatterer.vol_frac = 0.001; op = VtsSpectroscopy.GetOP(absorbers, scatterer, wv); @@ -384,7 +383,8 @@ xlabel('Wavelength, \lambda [nm]'); options = [{'Location', 'NorthWest'}; {'FontSize', 12}; {'Box', 'on'}]; PlotHelper.CreateLegend(fx, 'f_x = ', 'mm^-^1', options); -%% Example ROfFx (multiple optical properties, varying mua as a function of wavelength, mus' as mie scatterer) + +%% Example 15: ROfFx (multiple optical properties, varying mua as a function of wavelength, mus' as mie scatterer) % Evaluate reflectance as a function of spacial-frequency with multiple % sets of optical properties, varying mua as a function of wavelength. @@ -403,14 +403,14 @@ % % or % scatterer.Type = 'Intralipid'; -% scatterer.vol_frac = 0.5; +% scatterer.vol_frac = 0.01; % % or scatterer.Type = 'Mie'; scatterer.radius = 0.5; scatterer.n = 1.4; scatterer.nMedium = 1.0; -scatterer.vol_frac = 0.5; +scatterer.vol_frac = 0.001; op = VtsSpectroscopy.GetOP(absorbers, scatterer, wv); @@ -445,7 +445,7 @@ options = [{'Location', 'NorthWest'}; {'FontSize', 12}; {'Box', 'on'}]; PlotHelper.CreateLegend(fx, 'f_x = ', 'mm^-^1', options); -%% Example ROfFx +%% Example 16: ROfFx % Call planar reflectance with multiple sets of optical % properties, varying the scattering prefactor as a function of wavelength. @@ -478,7 +478,7 @@ options = [{'Location', 'NorthWest'}; {'FontSize', 12}; {'Box', 'on'}]; PlotHelper.CreateLegend(A, '\mu_s''(1000nm) = ', 'mm^-^1', options); -%% Example ROfRho (multiple wavelengths, multiple rho) +%% Example 17: ROfRho (multiple wavelengths, multiple rho) % Call reflectance varying the wavelength. VtsSolvers.SetSolverType('PointSourceSDA'); @@ -508,7 +508,7 @@ xlabel('Wavelength, \lambda [nm]'); options = [{'Location', 'NorthEast'}; {'FontSize', 12}; {'Box', 'on'}]; PlotHelper.CreateLegend(rho,'\rho = ', 'mm',options); -%% Example ROfRho (inverse solution for chromophore concentrations for multiple wavelengths, single rho) +%% Example 18: ROfRho (inverse solution for chromophore concentrations for multiple wavelengths, single rho) rho = 1; % source-detector separation in mm wv = 400:50:1000; @@ -571,7 +571,8 @@ (measData-recovered)*(measData-recovered)')); disp(sprintf('error = [%5.3f %5.3f %5.3f]',abs(measData(1)-recovered(1))/measData(1),... abs(measData(2)-recovered(2))/measData(2),abs(measData(3)-recovered(3))/measData(3))); -%% Example ROfRho for a two-layer tissue (multiple optical properties and rhos) + +%% Example 19: ROfRho for a two-layer tissue (multiple optical properties and rhos) clear op topLayerThickness = 2; % units: mm @@ -592,7 +593,8 @@ title('2-Layer Reflectance vs \rho for various top Layer OPs'); ylabel('R(\rho)'); xlabel('\rho'); -%% Example ROfFx for a two-layer tissue (multiple optical properties and fxs) + +%% Example 20: ROfFx for a two-layer tissue (multiple optical properties and fxs) clear op topLayerThickness = 2; % units: mm @@ -605,7 +607,7 @@ fx = 0:0.01:0.5; % range of spatial frequencies in 1/mm test = VtsSolvers.ROfFxTwoLayer(op, topLayerThickness, fx); -f = figure; semilogy(fx, test); +f = figure; plot(fx, test); set(f,'Name','2-Layer Reflectance vs Fx for various top Layer Optical Properties'); % create the legend with just the mua value from the top layer optical properties options = [{'FontSize', 12}; {'Location', 'NorthEast'}]; @@ -613,7 +615,8 @@ title('2-Layer Reflectance vs fx for various top Layer OPs'); ylabel('R(fx)'); xlabel('fx'); -%% Example ROfRhoAndTime for a two-layer tissue (multiple optical properties and times) + +%% Example 21: ROfRhoAndTime for a two-layer tissue (multiple optical properties and times) clear op topLayerThickness = 2; % units: mm % @@ -662,7 +665,7 @@ ylabel('R(t)'); xlabel('Time, t [ns]'); -%% Example ROfRhoAndFt for a two-layer tissue (multiple optical properties and fts) +%% Example 22: ROfRhoAndFt for a two-layer tissue (multiple optical properties and fts) % Evaluate reflectance as a function of rho and temporal-frequency with one % set of optical properites. clear op diff --git a/src/Vts.Scripting.Test/GlobalUsings.cs b/src/Vts.Scripting.Test/GlobalUsings.cs new file mode 100644 index 000000000..cefced496 --- /dev/null +++ b/src/Vts.Scripting.Test/GlobalUsings.cs @@ -0,0 +1 @@ +global using NUnit.Framework; \ No newline at end of file diff --git a/src/Vts.Scripting.Test/ProgramTests.cs b/src/Vts.Scripting.Test/ProgramTests.cs new file mode 100644 index 000000000..2ef8248b4 --- /dev/null +++ b/src/Vts.Scripting.Test/ProgramTests.cs @@ -0,0 +1,27 @@ +namespace Vts.Scripting.Test; + +public class ProgramTests +{ + [SetUp] + public void Setup() + { + } + + [Test] + public void Confirm_RunAllMonteCarloDemos_Does_Not_Throw() + { + Assert.DoesNotThrow(() => BatchDemoRunner.RunAllMonteCarloDemos(showPlots: false)); + } + + [Test] + public void Confirm_RunAllForwardSolverDemos_Does_Not_Throw() + { + Assert.DoesNotThrow(() => BatchDemoRunner.RunAllForwardSolverDemos(showPlots: false)); + } + + [Test] + public void Confirm_RunAllShortCourseDemos_Does_Not_Throw() + { + Assert.DoesNotThrow(() => BatchDemoRunner.RunAllShortCourseDemos(showPlots: false)); + } +} \ No newline at end of file diff --git a/src/Vts.Scripting.Test/Vts.Scripting.Test.csproj b/src/Vts.Scripting.Test/Vts.Scripting.Test.csproj new file mode 100644 index 000000000..0fffbec3a --- /dev/null +++ b/src/Vts.Scripting.Test/Vts.Scripting.Test.csproj @@ -0,0 +1,26 @@ + + + + net6.0 + enable + enable + true + 11 + false + true + + + + + + + + + + + + + + + + diff --git a/src/Vts.Scripting/BatchDemoRunner.cs b/src/Vts.Scripting/BatchDemoRunner.cs new file mode 100644 index 000000000..4072ba441 --- /dev/null +++ b/src/Vts.Scripting/BatchDemoRunner.cs @@ -0,0 +1,64 @@ +[assembly: System.Runtime.CompilerServices.InternalsVisibleTo("Vts.Scripting.Test")] +namespace Vts.Scripting; + +/// +/// Static class to run all demos +/// +internal static class BatchDemoRunner +{ + /// + /// Method to run all Monte Carlo demos + /// + internal static void RunAllMonteCarloDemos(bool showPlots = true) + { + MonteCarlo.Demo01ROfRhoSimple.RunDemo(showPlots); + MonteCarlo.Demo02DawVsCaw.RunDemo(showPlots); + MonteCarlo.Demo03ROfRhoFullCustomization.RunDemo(showPlots); + MonteCarlo.Demo04N1000VsN100.RunDemo(showPlots); + MonteCarlo.Demo05PostProcessor.RunDemo(showPlots); + MonteCarlo.Demo06PmcPostProcessor.RunDemo(showPlots); + MonteCarlo.Demo07PmcInversion.RunDemo(showPlots); + MonteCarlo.Demo08UnitTestComparison.RunDemo(showPlots); + MonteCarlo.Demo09TransmittanceTallies.RunDemo(showPlots); + MonteCarlo.Demo10ROfFx.RunDemo(showPlots); + } + + /// + /// Method to run all Forward Solver demos + /// + internal static void RunAllForwardSolverDemos(bool showPlots = true) + { + ForwardSolvers.Demo01ROfRhoAndFtSingle.RunDemo(showPlots); + ForwardSolvers.Demo02ROfFxAndFtMulti.RunDemo(showPlots); + ForwardSolvers.Demo03FluenceOfRhoAndZ.RunDemo(showPlots); + ForwardSolvers.Demo04FluenceOfRhoAndZAndFt.RunDemo(showPlots); + ForwardSolvers.Demo05PhdOfRhoAndZ.RunDemo(showPlots); + ForwardSolvers.Demo06FluenceOfRhoAndZTwoLayer.RunDemo(showPlots); + ForwardSolvers.Demo07PhdOfRhoAndZTwoLayer.RunDemo(showPlots); + ForwardSolvers.Demo08AbsorbedEnergyOfRhoAndZ.RunDemo(showPlots); + ForwardSolvers.Demo09ROfRhoMulti.RunDemo(showPlots); + ForwardSolvers.Demo10ROfRhoAndTimeMulti.RunDemo(showPlots); + ForwardSolvers.Demo11ROfFxAndTime.RunDemo(showPlots); + ForwardSolvers.Demo12ROfFxSingle.RunDemo(showPlots); + ForwardSolvers.Demo13ROfFxMultiOpProp.RunDemo(showPlots); + ForwardSolvers.Demo14ROfFxMultiOpPropIntralipid.RunDemo(showPlots); + ForwardSolvers.Demo15ROfFxMultiOpPropMie.RunDemo(showPlots); + ForwardSolvers.Demo16ROfFxMultiPowerLaw.RunDemo(showPlots); + ForwardSolvers.Demo17ROfRhoMultiOpProp.RunDemo(showPlots); + ForwardSolvers.Demo18ROfRhoMultiOpPropInversion.RunDemo(showPlots); + ForwardSolvers.Demo19ROfRhoTwoLayerMultiOpProp.RunDemo(showPlots); + ForwardSolvers.Demo20ROfFxTwoLayerMultiOpProp.RunDemo(showPlots); + ForwardSolvers.Demo21ROfRhoAndFtTwoLayerMultiOpProp.RunDemo(showPlots); + ForwardSolvers.Demo22ROfRhoAndFtTwoLayerMultiOpProp.RunDemo(showPlots); + } + + /// + /// Method to run all Short Course examples + /// + internal static void RunAllShortCourseDemos(bool showPlots = true) + { + ShortCourse.Demo01APhotonCountWithFluence.RunDemo(showPlots); + ShortCourse.Demo01BAnalogVsDiscreteWithFluence.RunDemo(showPlots); + ShortCourse.Demo02AnalogVsContinuousWithReflectance.RunDemo(showPlots); + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo01ROfRhoAndFtSingle.cs b/src/Vts.Scripting/ForwardSolvers/Demo01ROfRhoAndFtSingle.cs new file mode 100644 index 000000000..a64824df5 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo01ROfRhoAndFtSingle.cs @@ -0,0 +1,44 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of rho and time frequency +/// +internal class Demo01ROfRhoAndFtSingle : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 01: predict R(rho, ft) based on a standard diffusion approximation solution to the time-dependent RTE + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var op = new OpticalProperties(mua: 0.01, musp: 1.2, g: 0.8, n: 1.4); + var rho = 10; // s-d separation, in mm + var fts = new DoubleRange(start: 0, stop: 0.5, number: 51).AsEnumerable().ToArray(); // range of temporal frequencies in GHz + + // predict the temporal frequency response at the specified optical properties and s-d separation + var rOfFt = solver.ROfRhoAndFt(op, rho, fts); + var magnitudeScale = 1E4; + var rOfFtAmplitude = rOfFt.Select(r => r.Magnitude * magnitudeScale).ToArray(); + var rOfFtPhase = rOfFt.Select(r => -r.Phase).ToArray(); + + //%% Example ROfRhoAndFt + //% Evaluate reflectance as a function of rho and temporal-frequency with one + var xLabel = "time frequency [GHz]"; + var charts = new[] + { + LineChart(fts, rOfFtAmplitude, xLabel, yLabel: $"|R(ft)@ρ={rho}mm| [mm-2]*{magnitudeScale:E0}"), + LineChart(fts, rOfFtPhase, xLabel, yLabel: $"Φ(R(ft)@ρ={rho}mm) [rad]") + }; + var grid = Chart.Grid(charts, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + if(showPlots) + { + grid.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo02ROfFxAndFtMulti.cs b/src/Vts.Scripting/ForwardSolvers/Demo02ROfFxAndFtMulti.cs new file mode 100644 index 000000000..5144e0b88 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo02ROfFxAndFtMulti.cs @@ -0,0 +1,56 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency and time frequency +/// +internal class Demo02ROfFxAndFtMulti : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 02: predict R(fx, ft) based on a standard diffusion approximation solution to the time-dependent RTE + // Note: + // - this example computes and displays R(ft) at multiple spatial frequencies, to demonstrate that multiplexing functionality + // - updating the fx range to emit a single value will isolate to the frequency desired + // - e.g. var fxs = new DoubleRange(start: 0, stop: 0, number: 1) will plot R(ft) at fx=0 only (i.e. planar illumination) + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var op = new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4); + var fxs = new DoubleRange(start: 0, stop: 0.1, number: 4).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + var fts = new DoubleRange(start: 0, stop: 0.5, number: 51).AsEnumerable().ToArray(); // range of temporal frequencies in GHz + + // predict the temporal frequency response at each specified optical property and spatial frequency + var rOfFxAndFt = fxs.Select(fx => solver.ROfFxAndFt(op, fx, fts)).ToArray(); + + // Plot reflectance as a function of temporal-frequency at the specified spatial frequencies + var xLabel = "time frequency [GHz]"; + var amplitudeAndPhaseCharts = new[] + { + Chart.Combine(fxs.Select((fx, fxi) => + LineChart(fts, rOfFxAndFt[fxi].Select(r => r.Magnitude).ToArray(), xLabel, yLabel: $"|R(ft)| [unitless]", title: $"amp@fx={fx:F3}"))), + Chart.Combine(fxs.Select((fx, fxi) => + LineChart(fts, rOfFxAndFt[fxi].Select(r => -r.Phase).ToArray(), xLabel, yLabel: $"Φ(R(ft)) [rad]", title: $"phase@fx={fx:F3}"))) + }; + var grid1 = Chart.Grid(amplitudeAndPhaseCharts, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + var realAndImaginaryCharts = new[] + { + Chart.Combine(fxs.Select((fx, fxi) => + LineChart(fts, rOfFxAndFt[fxi].Select(r => r.Real).ToArray(), xLabel, yLabel: $"R(ft) (real) [unitless]", title: $"real@fx={fx:F3}"))), + Chart.Combine(fxs.Select((fx, fxi) => + LineChart(fts, rOfFxAndFt[fxi].Select(r => r.Imaginary).ToArray(), xLabel, yLabel: $"R(ft) (imag) [unitless]", title: $"imag@fx={fx:F3}"))) + }; + var grid2 = Chart.Grid(realAndImaginaryCharts, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + if (showPlots) + { + grid1.Show(); + grid2.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo03FluenceOfRhoAndZ.cs b/src/Vts.Scripting/ForwardSolvers/Demo03FluenceOfRhoAndZ.cs new file mode 100644 index 000000000..f4afafb5b --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo03FluenceOfRhoAndZ.cs @@ -0,0 +1,66 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of +/// radial extent and depth using optical properties from a list of chromophore absorbers +/// with their concentrations and a power law scatterer for a range of wavelengths +/// +internal class Demo03FluenceOfRhoAndZ : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 03: Evaluate fluence as a function of radial extent and depth using optical properties from a list of + // chromophore absorbers with their concentrations and a power law scatterer for a range of wavelengths. + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var rhos = new DoubleRange(start: 0.1, stop: 10, number: 100).AsEnumerable().ToArray(); // range of s-d separations in mm + var zs = new DoubleRange(start: 0.1, stop: 10, number: 100).AsEnumerable().ToArray(); // range of depths in mm + + // retrieve desired optical properties, based on spectral data information + + // create an array of chromophore absorbers, each with a given concentrations + var chromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, 70), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, 30), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, 0.8), // fractional concentration + }; + + // construct a scatterer + var scatterer = new PowerLawScatterer(a: 1.2, b: 1.42); + // or: var scatterer = new IntralipidScatterer(volumeFraction: 0.5); + // or: var scatterer = new MieScatterer(particleRadius: 0.5, particleRefractiveIndex: 1.4, mediumRefractiveIndex: 1.0, volumeFraction: 0.5); + + // compose a tissue using the chromophores and scatterer + var tissue = new Tissue(chromophores, scatterer, "", n: 1.4); + + // predict the tissue's fluence(rho, z) for tissue optical properties spanning the visible and NIR spectral regimes + var wavelengths = new DoubleRange(start: 450, stop: 1000, number: 1101).AsEnumerable().ToArray(); // range of wavelengths in nm + var op = tissue.GetOpticalProperties(wavelengths); + var fluenceOfRhoAndZ = solver.FluenceOfRhoAndZ(op, rhos, zs); + + // Plot the log(fluence(rho, z)) for the last wavelength + var wvi = wavelengths.Length - 1; // index of last wavelength (1000nm) + var imageSize = rhos.Length * zs.Length; + var allRhos = rhos.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + var fluenceRowsToPlot = fluenceOfRhoAndZ + .Skip(wvi * imageSize) // skip to the last wavelength + .Select(fluence => Math.Log(fluence)) // take log for visualization purposes + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var fluenceDataToPlot = fluenceRowsToPlot.Reverse().Concat(fluenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var map = Heatmap(values: fluenceDataToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"log(Φ(ρ, z) @λ={wavelengths[wvi]}nm"); + + if (showPlots) + { + map.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo04FluenceOfRhoAndZAndFt.cs b/src/Vts.Scripting/ForwardSolvers/Demo04FluenceOfRhoAndZAndFt.cs new file mode 100644 index 000000000..48cfb18d9 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo04FluenceOfRhoAndZAndFt.cs @@ -0,0 +1,70 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting fluence as a function of +/// radial extent, depth, and time frequency at a given set of optical properties +/// +internal class Demo04FluenceOfRhoAndZAndFt : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 04: Evaluate fluence as a function of radial extent, depth, + // and time frequency at a given set of optical properties + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var op = new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4); + var rhos = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of s-d separations in mm + var zs = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of depths in mm + var fts = new DoubleRange(start: 0, stop: 1, number: 2).AsEnumerable().ToArray(); // range of time frequencies in GHz + + // predict the tissue's fluence(rho, z, ft) for the given optical properties + var fluenceOfRhoAndZAndFt = solver.FluenceOfRhoAndZAndFt(new[]{ op }, rhos, zs, fts); + + // collect the steady-state (CW) and frequency domain fluence maps for plotting + // (innermost axis is ft, so we need to take every 2nd starting at index 0) + var (ft1, ft2) = (0, 1); // indices of the two time frequencies, ft, to plot + var cwFluenceOfRhoAndZ = fluenceOfRhoAndZAndFt.TakeEveryNth(n: fts.Length, skip: ft1).ToArray(); + var fdFluenceOfRhoAndZ = fluenceOfRhoAndZAndFt.TakeEveryNth(n: fts.Length, skip: ft2).ToArray(); + + var allRhos = rhos.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + // Plot the CW fluence: log(fluence(rho, z, ft=0GHz)) + var cwFluenceRowsToPlot = cwFluenceOfRhoAndZ + .Select(fluence => Math.Log(fluence.Magnitude)) // take log for visualization purposes + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var allCwFluenceRowsToPlot = cwFluenceRowsToPlot.Reverse().Concat(cwFluenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var cwFluenceChart = Heatmap(values: allCwFluenceRowsToPlot, x: allRhos, y: zs, xLabel: "ρ", yLabel: "z", title: "log(Φ(ρ, z, ft=0Ghz))"); + + // Plot the frequency-domain fluence amplitude at 1GHz: log(fluence(rho, z, ft=0GHz)) + var fdFluenceRowsToPlot = fdFluenceOfRhoAndZ + .Select(fluence => Math.Log(fluence.Magnitude)) // take log for visualization purposes + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var allFdFluenceRowsToPlot = fdFluenceRowsToPlot.Reverse().Concat(fdFluenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var fdFluenceChart = Heatmap(values: allFdFluenceRowsToPlot, x: allRhos, y: zs, xLabel: "ρ", yLabel: "z", title: "log(Φ(ρ, z, ft=1Ghz))"); + + // calculate the modulation by taking an element-wise ratio of the fluence maps (i.e. FD/Cw) + var modFluenceOfRhoAndZ = fdFluenceOfRhoAndZ.Zip(cwFluenceOfRhoAndZ, (fd, cw) => fd.Magnitude / cw.Magnitude).ToArray(); + + // Plot the modulation at 1GHz: log(fluence(rho, z, ft=1GHz)/fluence(rho, z, ft=0GHz)) + var modFluenceRowsToPlot = modFluenceOfRhoAndZ + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var allModFluenceRowsToPlot = modFluenceRowsToPlot.Reverse().Concat(modFluenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var modFluenceChart = Heatmap(values: allModFluenceRowsToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: "modulation(ρ, z) @ ft=1Ghz))"); + + if (showPlots) + { + cwFluenceChart.Show(); + fdFluenceChart.Show(); + modFluenceChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo05PHDOfRhoAndZ.cs b/src/Vts.Scripting/ForwardSolvers/Demo05PHDOfRhoAndZ.cs new file mode 100644 index 000000000..6e9d7718c --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo05PHDOfRhoAndZ.cs @@ -0,0 +1,47 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting photon hitting density (PHD) +/// as a function of radial extent and depth at a given set of optical properties and source-detector separation +/// +internal class Demo05PhdOfRhoAndZ : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 05: Compute photon hitting density (PHD) as a function of radial extent + // and depth at a given set of optical properties and source-detector separation + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var op = new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4); + var rhos = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of s-d separations in mm + var zs = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of depths in mm + + // predict the tissue's fluence(rho, z) for the given optical properties + var fluenceOfRhoAndZ = solver.FluenceOfRhoAndZ(new[]{ op }, rhos, zs); + var sourceDetectorSeparation = 10; // mm + var phd = ComputationFactory.GetPHD(forwardSolverType: ForwardSolverType.PointSourceSDA, + fluenceOfRhoAndZ, sourceDetectorSeparation, ops: new[] { op }, rhos, zs); + + var allRhos = rhos;//.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + // Plot the CW fluence: log(fluence(rho, z, ft=0GHz)) + var phdRowsToPlot = phd + .Select(fluence => Math.Log(fluence)) // take log for visualization purposes + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var allPhdRowsToPlot = phdRowsToPlot;//.Reverse().Concat(phdRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var phdChart = Heatmap( + values: allPhdRowsToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"PHD(ρ, z) @ s-d: {sourceDetectorSeparation} mm"); + + if (showPlots) + { + phdChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo06FluenceOfRhoAndZTwoLayer.cs b/src/Vts.Scripting/ForwardSolvers/Demo06FluenceOfRhoAndZTwoLayer.cs new file mode 100644 index 000000000..7e897f28d --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo06FluenceOfRhoAndZTwoLayer.cs @@ -0,0 +1,56 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting fluence in a multi-layer tissue +/// as a function of radial extent and depth at a given set of optical properties +/// +internal class Demo06FluenceOfRhoAndZTwoLayer : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 06: Compute fluence for a two-layer medium as a function + // of radial extent and depth at a given set of optical properties + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new TwoLayerSDAForwardSolver { SourceConfiguration = SourceConfiguration.Distributed }; + var topLayerThickness = 5; // mm + var opRegions = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion + ( + zRange: new DoubleRange(0, topLayerThickness, 2), + regionOP: new OpticalProperties(mua: 0.1, musp: 1, g: 0.8, n: 1.4) + ), + new LayerOpticalPropertyRegion + ( + zRange: new DoubleRange(topLayerThickness, double.PositiveInfinity, 2), + regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4) + ) + }; + var rhos = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of s-d separations in mm + var zs = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of depths in mm + + // predict the tissue's fluence(rho, z) for the given optical properties + var fluenceOfRhoAndZ = solver.FluenceOfRhoAndZ(new[] { opRegions }, rhos, zs ); + + var allRhos = rhos.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + var fluenceRowsToPlot = fluenceOfRhoAndZ + .Select(fluence => Math.Log(fluence)) // take log for visualization purposes + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var allFluenceRowsToPlot = fluenceRowsToPlot.Reverse().Concat(fluenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var fluenceChart = Heatmap( + values: allFluenceRowsToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"log(Φ(ρ, z) [mm-3])"); + + if (showPlots) + { + fluenceChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo07PHDOfRhoAndZTwoLayer.cs b/src/Vts.Scripting/ForwardSolvers/Demo07PHDOfRhoAndZTwoLayer.cs new file mode 100644 index 000000000..175670962 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo07PHDOfRhoAndZTwoLayer.cs @@ -0,0 +1,61 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting photon hitting density (PHD) for a two-layer medium +/// as a function of radial extent and depth at a given set of optical properties and source-detector separation +/// +internal class Demo07PhdOfRhoAndZTwoLayer : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 07: Compute photon hitting density (PHD) for a two-layer medium as a function + // of radial extent and depth at a given set of optical properties and source-detector separation + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new TwoLayerSDAForwardSolver { SourceConfiguration = SourceConfiguration.Distributed }; + var topLayerThickness = 5; // mm + var opRegions = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion + ( + zRange: new DoubleRange(0, topLayerThickness, 2), + regionOP: new OpticalProperties(mua: 0.1, musp: 1, g: 0.8, n: 1.4) + ), + new LayerOpticalPropertyRegion + ( + zRange: new DoubleRange(topLayerThickness, double.PositiveInfinity, 2), + regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4) + ) + }; + var rhos = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of s-d separations in mm + var zs = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of depths in mm + + // predict the tissue's fluence(rho, z) for the given optical properties + var fluenceOfRhoAndZ = solver.FluenceOfRhoAndZ(new[] { opRegions }, rhos, zs ); + + var sourceDetectorSeparation = 10; // mm + var phd = ComputationFactory.GetPHD(forwardSolverType: ForwardSolverType.PointSourceSDA, + fluenceOfRhoAndZ, sourceDetectorSeparation, new[]{ opRegions[0].RegionOP }, rhos, zs); // todo: pick op based on layer thickness? + + // plot the PHD @ 10 mm s-d separation + var allRhos = rhos;//.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + var phdRowsToPlot = phd + .Select(fluence => Math.Log(fluence)) // take log for visualization purposes + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var allPhdRowsToPlot = phdRowsToPlot;//.Reverse().Concat(phdRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var phdChart = Heatmap( + values: allPhdRowsToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"PHD(ρ, z) [unitless] @ s-d: {sourceDetectorSeparation} mm"); + + if (showPlots) + { + phdChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo08AbsorbedEnergyOfRhoAndZ.cs b/src/Vts.Scripting/ForwardSolvers/Demo08AbsorbedEnergyOfRhoAndZ.cs new file mode 100644 index 000000000..e023ea60c --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo08AbsorbedEnergyOfRhoAndZ.cs @@ -0,0 +1,44 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting absorbed energy as a function of +/// radial extent and depth at a given set of optical properties +/// +internal class Demo08AbsorbedEnergyOfRhoAndZ : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 08: Compute absorbed energy as a function of radial extent and depth at a given set of optical properties + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver { SourceConfiguration = SourceConfiguration.Distributed }; + var op = new OpticalProperties(mua: 0.1, musp: 1, g: 0.8, n: 1.4); + var rhos = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of s-d separations in mm + var zs = new DoubleRange(start: 0.1, stop: 19.9, number: 100).AsEnumerable().ToArray(); // range of depths in mm + + // predict the tissue's AbsorbedEnergy(rho, z) for the given optical properties + var fluenceOfRhoAndZ = solver.FluenceOfRhoAndZ(new[] { op }, rhos, zs ); + var absorbedEnergyOfRhoAndZ = ComputationFactory.GetAbsorbedEnergy(fluenceOfRhoAndZ, op.Mua).ToArray(); + + var allRhos = rhos.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + var absorbedEnergyRowsToPlot = absorbedEnergyOfRhoAndZ + .Select(ae => Math.Log(ae)) // take log for visualization purposes + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var allAbsorbedEnergyRowsToPlot = absorbedEnergyRowsToPlot.Reverse() + .Concat(absorbedEnergyRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var absorbedEnergyChart = Heatmap( + values: allAbsorbedEnergyRowsToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"log(AbsorbedEnergy(ρ, z) [mm-3])"); + + if (showPlots) + { + absorbedEnergyChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo09ROfRhoMulti.cs b/src/Vts.Scripting/ForwardSolvers/Demo09ROfRhoMulti.cs new file mode 100644 index 000000000..b233e8786 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo09ROfRhoMulti.cs @@ -0,0 +1,42 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of source-detector separation +/// for multiple sets of optical properties +/// +internal class Demo09ROfRhoMulti : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 08: predict R(rho) for multiple OPs based on a standard diffusion approximation solution to the time-dependent RTE + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var ops = new[] + { + new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4), + new OpticalProperties(mua: 0.1, musp: 1, g: 0.8, n: 1.4), + new OpticalProperties(mua: 1, musp: 1, g: 0.8, n: 1.4) + }; + var rhos = new DoubleRange(start: 0.5, stop: 9.5, number: 19).AsEnumerable().ToArray(); // range of radial detector locations in mm + + // predict the reflectance at each specified optical property and source-detector separation + var allROfRho = solver.ROfRho(ops, rhos).ToArray(); + + // Plot log(reflectance) as a function of radial distance at the specified spatial frequencies + var allCharts = allROfRho.Chunk(rhos.Length).Select((rOfRho, rIdx) => // pull out each R(ρ) (outer loop is optical properties) + LineChart(rhos, rOfRho.Select(r => Math.Log(r)).ToArray(), + xLabel: "ρ [mm]", yLabel: $"log(R(ρ) [mm-2])", title: $"log(R(ρ)) @ mua={ops[rIdx].Mua:F3}")); + var chartCombined = Chart.Combine(allCharts); + + if (showPlots) + { + chartCombined.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo10ROfRhoAndTimeMulti.cs b/src/Vts.Scripting/ForwardSolvers/Demo10ROfRhoAndTimeMulti.cs new file mode 100644 index 000000000..cbe84b49d --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo10ROfRhoAndTimeMulti.cs @@ -0,0 +1,44 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of source-detector separation +/// and time for multiple sets of optical properties +/// +internal class Demo10ROfRhoAndTimeMulti : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 10: predict R(rho,t) at various OPs based on a standard diffusion approximation solution to the time-dependent RTE + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var ops = new[] + { + new OpticalProperties(mua: 0.1, musp: 1, g: 0.8, n: 1.4), + new OpticalProperties(mua: 0.2, musp: 1, g: 0.8, n: 1.4) + }; + var rho = 10; // radial detector location in mm + var ts = new DoubleRange(start: 0, stop: 0.5, number: 501).AsEnumerable().ToArray(); // range of times in ns + + // predict the reflectance at each specified optical property and source-detector separation + var allROfRho = solver.ROfRhoAndTime(ops, rho, ts).ToArray(); + + // Plot log(reflectance) as a function of time at a range of source-detector separations + // Create two plots - one for each set of optical properties + var bothCharts = allROfRho.Chunk(ts.Length).Select((rOfTime, opIdx) => // pull out each R(t) (innermost loop is time) + LineChart(ts, rOfTime.Select(r =>r).ToArray(), + xLabel: "t [ns]", yLabel: $"R(t) [mm-2*ns-1])", title: $"R(t) @ mua={ops[opIdx].Mua:F3}")); + + var chartCombined = Chart.Combine(bothCharts); + + if (showPlots) + { + chartCombined.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo11ROfFxAndTime.cs b/src/Vts.Scripting/ForwardSolvers/Demo11ROfFxAndTime.cs new file mode 100644 index 000000000..77f51415d --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo11ROfFxAndTime.cs @@ -0,0 +1,47 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency and time +/// +internal class Demo11ROfFxAndTime : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 11: predict R(fx, t) based on a standard diffusion approximation solution to the time-dependent RTE + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var op = new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4); + var fxs = new DoubleRange(start: 0, stop: 0.5, number: 11).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + var ts = new DoubleRange(start: 0, stop: 0.5, number: 501).AsEnumerable().ToArray(); // range of times in ns + + // predict the temporal response at each specified optical property and spatial frequency + var rOfFxAndTime = solver.ROfFxAndTime(op, fxs, ts).ToArray(); + + // Plot reflectance as a function of time at the specified spatial frequencies + var fxChart = Chart.Combine(rOfFxAndTime.Chunk(ts.Length).Select((rOfTime ,fxi) => // take and plot each R(t) (outer loop is spatial frequency) + LineChart(ts, rOfTime, xLabel: "time [ns]", yLabel: $"R(t) [ns-1]", title: $"amp@fx={fxs[fxi]:F3}"))); + + var (tIdx0p01, tIdx0p05) = (10, 50); + var rOfFx0p01 = rOfFxAndTime.TakeEveryNth(ts.Length, skip: tIdx0p01).ToArray(); // t = 0.01 ns data + var rOfFx0p05 = rOfFxAndTime.TakeEveryNth(ts.Length, skip: tIdx0p05).ToArray(); // t = 0.05 ns data + + // Plot reflectance as a function of time at the specified spatial frequencies + var timeChart = Chart.Combine(new[] + { + LineChart(fxs, rOfFx0p01, xLabel: "fx [mm-1]", yLabel: $"R(fx) [unitless]", title: $"amp@t={ts[tIdx0p01]:F3}"), + LineChart(fxs, rOfFx0p05, xLabel: "fx [mm-1]", yLabel: $"R(fx) [unitless]", title: $"amp@t={ts[tIdx0p05]:F3}") + }); + + if (showPlots) + { + fxChart.Show(); + timeChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo12ROfFxSingle.cs b/src/Vts.Scripting/ForwardSolvers/Demo12ROfFxSingle.cs new file mode 100644 index 000000000..22f884900 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo12ROfFxSingle.cs @@ -0,0 +1,34 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency at a single set of optical properties +/// +internal class Demo12ROfFxSingle : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 12: predict R(fx) based on a standard diffusion approximation solution to the time-dependent RTE + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var op = new OpticalProperties(mua: 0.1, musp: 1.2, g: 0.8, n: 1.4); + var fxs = new DoubleRange(start: 0, stop: 0.2, number: 201).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + + // predict the spatial frequency response at each specified optical properties + var rOfFx = solver.ROfFx(op, fxs); + + // Plot reflectance as a function of spatial frequency + var rOfFxAmplitude = rOfFx.Select(Math.Abs).ToArray(); + var chart = LineChart(fxs, rOfFxAmplitude, xLabel: "fx [mm-1]", yLabel: $"R(fx) [unitless]", title: "Reflectance vs spatial frequency"); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo13ROfFxMultiOpProp.cs b/src/Vts.Scripting/ForwardSolvers/Demo13ROfFxMultiOpProp.cs new file mode 100644 index 000000000..790246c9c --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo13ROfFxMultiOpProp.cs @@ -0,0 +1,38 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency while varying absorption linearly +/// +internal class Demo13ROfFxMultiOpProp : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 13: predict R(fx) based on a standard diffusion approximation solution to the time-dependent RTE + // while varying absorption linearly + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var fxs = new DoubleRange(start: 0, stop: 0.2, number: 201).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + var muas = new DoubleRange(start: 0, stop: 0.1, number: 11).AsEnumerable().ToArray(); // range of absorption values in 1/mm + var ops = muas.Select(mua => new OpticalProperties(mua: mua, musp: 1.2, g: 0.8, n: 1.4)).ToArray(); + + // predict the spatial frequency response at each specified optical properties + var rOfFx = solver.ROfFx(ops, fxs); + + // Plot reflectance as a function of spatial frequencies at each set of optical properties + var rOfFxAmplitude = rOfFx.Select(Math.Abs).ToArray(); + var chart = Chart.Combine(rOfFxAmplitude.Chunk(fxs.Length).Select((rOfFxSingleOpProp, opi) => // take and plot each R(fx) (outer loop is optical properties) + LineChart(fxs, rOfFxSingleOpProp, xLabel: "fx [mm-1]", yLabel: $"R(fx) [unitless] (varying mua linearly)", + title: $"Reflectance @ mua={ops[opi].Mua:F3}, musp={ops[opi].Musp:F3}"))); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo14ROfFxMultiOpPropIntralipid.cs b/src/Vts.Scripting/ForwardSolvers/Demo14ROfFxMultiOpPropIntralipid.cs new file mode 100644 index 000000000..f87f90bf1 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo14ROfFxMultiOpPropIntralipid.cs @@ -0,0 +1,65 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency +/// where optical properties vary as a function of wavelength using an intralipid scatterer +/// +internal class Demo14ROfFxMultiOpPropIntralipid : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 14: predict R(fx) based on a standard diffusion approximation solution to the time-dependent RTE + // where optical properties vary as a function of wavelength using an intralipid scatterer + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var fxs = new DoubleRange(start: 0, stop: 0.2, number: 5).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + + // retrieve desired optical properties, based on spectral data information + + // create an array of chromophore absorbers, each with a given concentrations + var chromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, 70), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, 30), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, 0.8), // fractional concentration + }; + + // construct a scatterer + var scatterer = new IntralipidScatterer(volumeFraction: 0.01); + // or: var scatterer = new MieScatterer(particleRadius: 0.5, particleRefractiveIndex: 1.4, mediumRefractiveIndex: 1.0, volumeFraction: 0.001); + // or: var scatterer = new PowerLawScatterer(a: 1.2, b: 1.42); + + // compose a tissue using the chromophores and scatterer + var tissue = new Tissue(chromophores, scatterer, "", n: 1.4); + + // predict the tissue's optical properties spanning the visible and NIR spectral regimes + var wavelengths = new DoubleRange(start: 450, stop: 1000, number: 1101).AsEnumerable().ToArray(); // range of wavelengths in nm + var ops = tissue.GetOpticalProperties(wavelengths); + + // predict the spatial frequency response at each specified optical properties + var rOfFx = solver.ROfFx(ops, fxs); + + // Plot mua, log(mua), musp, and reflectance as a function of wavelength at each set of optical properties + var rOfFxAmplitude = rOfFx.Select(Math.Abs).ToArray(); + var chart = Chart.Grid(new[] + { + LineChart(wavelengths, ops.Select(op => op.Mua).ToArray(), xLabel: "wavelength [nm]", yLabel: $"mua", title: "mua [mm-1]"), + LineChart(wavelengths, ops.Select(op => Math.Log(op.Mua)).ToArray(), xLabel: "wavelength [nm]", yLabel: $"log(mua)", title: "log(mua [mm-1])"), + LineChart(wavelengths, ops.Select(op => op.Musp).ToArray(), xLabel: "wavelength [nm]", yLabel: $"musp", title: "musp [mm-1]"), + Chart.Combine(fxs.Select((fx, fxi) => // plot R(wavelength)@fx for each spatial frequency, fx + LineChart(wavelengths, rOfFxAmplitude.TakeEveryNth(fxs.Length, skip: fxi).ToArray(), + xLabel: "wavelength [nm]", yLabel: $"R(wv)", title: $"R@fx={fx:F3} mm-1"))) + }, nRows: 4, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo15ROfFxMultiOpPropMie.cs b/src/Vts.Scripting/ForwardSolvers/Demo15ROfFxMultiOpPropMie.cs new file mode 100644 index 000000000..e2260624b --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo15ROfFxMultiOpPropMie.cs @@ -0,0 +1,65 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency +/// where optical properties vary as a function of wavelength using a Mie-theory scatterer +/// +internal class Demo15ROfFxMultiOpPropMie : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 15: predict R(fx) based on a standard diffusion approximation solution to the time-dependent RTE + // where optical properties vary as a function of wavelength using a Mie-theory scatterer + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var fxs = new DoubleRange(start: 0, stop: 0.2, number: 5).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + + // retrieve desired optical properties, based on spectral data information + + // create an array of chromophore absorbers, each with a given concentrations + var chromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, 70), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, 30), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, 0.8), // fractional concentration + }; + + // construct a scatterer + var scatterer = new MieScatterer(particleRadius: 0.5, particleRefractiveIndex: 1.4, mediumRefractiveIndex: 1.0, volumeFraction: 0.001); + // or: var scatterer = new IntralipidScatterer(volumeFraction: 0.01); + // or: var scatterer = new PowerLawScatterer(a: 1.2, b: 1.42); + + // compose a tissue using the chromophores and scatterer + var tissue = new Tissue(chromophores, scatterer, "", n: 1.4); + + // predict the tissue's optical properties spanning the visible and NIR spectral regimes + var wavelengths = new DoubleRange(start: 450, stop: 1000, number: 1101).AsEnumerable().ToArray(); // range of wavelengths in nm + var ops = tissue.GetOpticalProperties(wavelengths); + + // predict the spatial frequency response at each specified optical properties + var rOfFx = solver.ROfFx(ops, fxs); + + // Plot mua, log(mua), musp, and reflectance as a function of wavelength at each set of optical properties + var rOfFxAmplitude = rOfFx.Select(Math.Abs).ToArray(); + var chart = Chart.Grid(new[] + { + LineChart(wavelengths, ops.Select(op => op.Mua).ToArray(), xLabel: "wavelength [nm]", yLabel: $"mua", title: "mua [mm-1]"), + LineChart(wavelengths, ops.Select(op => Math.Log(op.Mua)).ToArray(), xLabel: "wavelength [nm]", yLabel: $"log(mua)", title: "log(mua [mm-1])"), + LineChart(wavelengths, ops.Select(op => op.Musp).ToArray(), xLabel: "wavelength [nm]", yLabel: $"musp", title: "musp [mm-1]"), + Chart.Combine(fxs.Select((fx, fxi) => // plot R(wavelength)@fx for each spatial frequency, fx + LineChart(wavelengths, rOfFxAmplitude.TakeEveryNth(fxs.Length, skip: fxi).ToArray(), + xLabel: "wavelength [nm]", yLabel: $"R(wv)", title: $"R@fx={fx:F3} mm-1"))) + }, nRows: 4, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo16ROfFxMultiPowerLaw.cs b/src/Vts.Scripting/ForwardSolvers/Demo16ROfFxMultiPowerLaw.cs new file mode 100644 index 000000000..c1b7e4aa1 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo16ROfFxMultiPowerLaw.cs @@ -0,0 +1,64 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency +/// where optical properties vary as a function of wavelength using a range of power-law scatterer prefactors, A +/// +internal class Demo16ROfFxMultiPowerLaw : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 16: predict R(fx) based on a standard diffusion approximation solution to the time-dependent RTE + // where optical properties vary as a function of wavelength using a range of power-law scatterer prefactors, A + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + + // retrieve desired optical properties, based on spectral data information + + // create an array of chromophore absorbers, each with a given concentrations + var chromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, 70), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, 30), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, 0.8), // fractional concentration + }; + + var fx = 0; // spatial frequency in 1/mm + + // predict the tissue's optical properties spanning the visible and NIR spectral regimes + var wavelengths = new DoubleRange(start: 450, stop: 1000, number: 1101).AsEnumerable().ToArray(); // range of wavelengths in nm + + var prefactorAs = new DoubleRange(start: 0.5, stop: 2.5, number: 9).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + var opsForMultipleA = prefactorAs.Select(prefactorA => + new Tissue(chromophores, new PowerLawScatterer(a: prefactorA, b: 1.42), "", n: 1.4).GetOpticalProperties(wavelengths)).ToArray(); + + // predict the wavelength response for each spectrum of optical properties given by a specific A prefactor + var rVsWavelengthForEachA = solver.ROfFx(opsForMultipleA.SelectMany(op => op).ToArray(), fx ); + + // Plot mua, log(mua), musp, and reflectance as a function of wavelength at each set of optical properties + var chart = Chart.Grid(new[] + { + LineChart(wavelengths, opsForMultipleA[0].Select(op => op.Mua).ToArray(), xLabel: "wavelength [nm]", yLabel: $"mua", title: "mua [mm-1]"), + LineChart(wavelengths, opsForMultipleA[0].Select(op => Math.Log(op.Mua)).ToArray(), xLabel: "wavelength [nm]", yLabel: $"log(mua)", title: "log(mua [mm-1])"), + Chart.Combine(prefactorAs.Select((prefactorA, pai) => + LineChart(wavelengths, opsForMultipleA[pai].Select(op => op.Musp).ToArray(), + xLabel: "wavelength [nm]", yLabel: $"musp", title: $"musp@A={prefactorA:F3} [mm-1]")) + ), + Chart.Combine(prefactorAs.Select((prefactorA, pai) => // plot R(wavelength) @ fx=0 for each prefactor A + LineChart(wavelengths, rVsWavelengthForEachA.Skip(wavelengths.Length * pai).Take(wavelengths.Length).ToArray(), + xLabel: "wavelength [nm]", yLabel: $"R(wv)", title: $"R@A={prefactorA:F3}")) + ) + }, nRows: 4, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo17ROfRhoMultiOpProp.cs b/src/Vts.Scripting/ForwardSolvers/Demo17ROfRhoMultiOpProp.cs new file mode 100644 index 000000000..7fd477b20 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo17ROfRhoMultiOpProp.cs @@ -0,0 +1,60 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of source-detector separation while varying absorption linearly +/// +internal class Demo17ROfRhoMultiOpProp : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 17: predict R(rho) based on a standard diffusion approximation solution to the time-dependent RTE + // while varying absorption linearly + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new PointSourceSDAForwardSolver(); + var rhos = new DoubleRange(start: 0.2, stop: 1, number: 5).AsEnumerable().ToArray(); // range of source-detector separations in mm + + // create an array of chromophore absorbers, each with a given concentrations + var chromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, 70), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, 30), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, 0.8), // fractional concentration + }; + + // construct a scatterer + var scatterer = new PowerLawScatterer(a: 1.2, b: 1.42); + + // compose a tissue using the chromophores and scatterer + var tissue = new Tissue(chromophores, scatterer, "", n: 1.4); + + // predict the tissue's optical properties spanning the visible and NIR spectral regimes + var wavelengths = new DoubleRange(start: 450, stop: 1000, number: 1101).AsEnumerable().ToArray(); // range of wavelengths in nm + var ops = tissue.GetOpticalProperties(wavelengths); + + // predict the radial reflectance response across the spectrum of optical properties + var rOfFx = solver.ROfRho(ops, rhos); + + // Plot mua, log(mua), musp, and reflectance as a function of wavelength at each set of optical properties + var rOfFxAmplitude = rOfFx.Select(Math.Abs).ToArray(); + var chart = Chart.Grid(new[] + { + LineChart(wavelengths, ops.Select(op => op.Mua).ToArray(), xLabel: "wavelength [nm]", yLabel: $"mua", title: "mua [mm-1]"), + LineChart(wavelengths, ops.Select(op => Math.Log(op.Mua)).ToArray(), xLabel: "wavelength [nm]", yLabel: $"log(mua)", title: "log(mua [mm-1])"), + LineChart(wavelengths, ops.Select(op => op.Musp).ToArray(), xLabel: "wavelength [nm]", yLabel: $"musp", title: "musp [mm-1]"), + Chart.Combine(rhos.Select((rho, rhoIdx) => // plot R(wavelength)@rho for each rho in mm + LineChart(wavelengths, rOfFxAmplitude.TakeEveryNth(rhos.Length, skip: rhoIdx).ToArray(), + xLabel: "wavelength [nm]", yLabel: $"R(wv)", title: $"R@rho={rho:F3} mm"))) + }, nRows: 4, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo18ROfRhoMultiOpPropInversion.cs b/src/Vts.Scripting/ForwardSolvers/Demo18ROfRhoMultiOpPropInversion.cs new file mode 100644 index 000000000..e0725b324 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo18ROfRhoMultiOpPropInversion.cs @@ -0,0 +1,114 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate using the Perturbation Monte Carlo post-processor to calculate optical properties (i.e. "inversion") +/// +internal class Demo18ROfRhoMultiOpPropInversion : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 18: determine chromophore concentration by performing a non-linear least squares fit to a simulated measured + // reflectance spectrum at a given source-detector separation, using a diffusion theory forward model + + // create an array of chromophore absorbers, each with a given concentrations + var chromophoresForMeasuredData = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, 70), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, 30), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, 0.8), // fractional concentration + }; + + // construct a scatterer + var scatterer = new PowerLawScatterer(a: 1.2, b: 1.42); + + // compose a tissue using the chromophores and scatterer + var tissue = new Tissue(chromophoresForMeasuredData, scatterer, "", n: 1.4); + + // predict the tissue's optical properties spanning the visible and NIR spectral regimes + //var wavelengths = new DoubleRange(start: 400, stop: 1000, number: 13).AsEnumerable().ToArray(); // range of wavelengths in nm + var wavelengths = new DoubleRange(start: 400, stop: 1000, number: 601).AsEnumerable().ToArray(); // range of wavelengths in nm + var measuredOPs = tissue.GetOpticalProperties(wavelengths); + + // Create some measurements, based on a Nurbs-based White Monte Carlo forward solver + var rho = 1; //source - detector separation in mm + var measurementForwardSolver = new NurbsForwardSolver(); + var measuredData = measurementForwardSolver.ROfRho(measuredOPs, rho); + + // Create a forward solver based on pMC prediction (see implementation below; note: implemented for ROfRho only) + var forwardSolverForInversion = new PointSourceSDAForwardSolver(); + + // declare a local forward reflectance function that computes reflectance from chromophores + // note that some variables are captured from the outer scope for simplicity (scatterer, wavelengths + double[] CalculateReflectanceVsWavelengthFromChromophoreConcentration(double[] chromophoreConcentration, params object[] otherValuesNeededForForwardSolution) + { + // create an array of chromophore absorbers, each with a given concentrations + var chromophoresLocal = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, chromophoreConcentration[0]), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, chromophoreConcentration[1]), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, chromophoreConcentration[2]), // fractional concentration + }; + + // compose a tissue using the chromophores and scatterer, and return its optical properties + var opsLocal = new Tissue(chromophoresLocal, scatterer, "", n: 1.4).GetOpticalProperties(wavelengths); + + // compute the reflectance based on that OP spectrum + var measuredDataLocal = forwardSolverForInversion.ROfRho(opsLocal, rho); + + return measuredDataLocal; + } + + // Run the inversion, based on Levenberg-Marquardt optimization. + // Note: chi-squared weighting in the inversion is based on standard deviation of measured data + var optimizer = new MPFitLevenbergMarquardtOptimizer(); + var initialGuess = new[] { 70, 30, 0.8 }; + var parametersToFit = new[] { true, true, true }; // fit all three chromophores simultaneously + var fit = optimizer.Solve( + initialGuess.Select(ig => ig).ToArray(), // make a copy, because the solver mutates this array + parametersToFit, + measuredData, + measuredData.Select(_ => 1.0).ToArray(), // weight spectrum equally + CalculateReflectanceVsWavelengthFromChromophoreConcentration); + + // calculate the final reflectance prediction at the fit values + var fitChromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, fit[0]), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, fit[1]), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, fit[2]), // fractional concentration + }; + var fitTissue = new Tissue(fitChromophores, scatterer, "", n: 1.4); + var fitOpticalProperties = fitTissue.GetOpticalProperties(wavelengths); + var fitReflectanceSpectrum = forwardSolverForInversion.ROfRho(fitOpticalProperties, rho); + + // calculate the initial guess reflectance prediction + var initialGuessChromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, initialGuess[0]), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, initialGuess[1]), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, initialGuess[2]), // fractional concentration + }; + var initialGuessTissueOpProp = new Tissue(initialGuessChromophores, scatterer, "", n: 1.4).GetOpticalProperties(wavelengths); + var initialGuessReflectanceSpectrum = forwardSolverForInversion.ROfRho(initialGuessTissueOpProp, rho); + + // plot and compare the results using Plotly.NET + var logMeasuredReflectance = measuredData.Select(r => Math.Log(r)).ToArray(); + var logGuessReflectance = initialGuessReflectanceSpectrum.Select(r => Math.Log(r)).ToArray(); + var logFitReflectance = fitReflectanceSpectrum.Select(r => Math.Log(r)).ToArray(); + var (xLabel, yLabel) = ("wavelength [nm]", "log(R(λ)) [mm-2]"); + var chart = Chart.Combine(new[] + { + ScatterChart(wavelengths, logMeasuredReflectance, xLabel, yLabel, title: "Measured Data"), + LineChart(wavelengths, logGuessReflectance, xLabel, yLabel, title: $"Initial guess (HbO2:{initialGuess[0]:F3}uM, Hb:{initialGuess[1]:F3}uM, H2O:{initialGuess[2]*100:F3}%)"), + LineChart(wavelengths, logFitReflectance, xLabel, yLabel, title: $"Converged result (HbO2:{fit[0]:F3}uM, Hb:{fit[1]:F3}uM, H2O:{fit[2]*100:F3}%)") + }); // show all three charts together + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo19ROfRhoTwoLayerMultiOpProp.cs b/src/Vts.Scripting/ForwardSolvers/Demo19ROfRhoTwoLayerMultiOpProp.cs new file mode 100644 index 000000000..dab059370 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo19ROfRhoTwoLayerMultiOpProp.cs @@ -0,0 +1,60 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of source-detector separation +/// for multiple sets of optical properties using a two-layer forward solver +/// +internal class Demo19ROfRhoTwoLayerMultiOpProp : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 19: predict R(rho) based on a standard diffusion approximation solution to the time-dependent RTE + // for multiple sets of optical properties using a two-layer forward solver + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new TwoLayerSDAForwardSolver(); + var rhos = new DoubleRange(start: 0.5, stop: 9.5, number: 19).AsEnumerable().ToArray(); // range of radial distances in 1/mm + var op1 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + var op2 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.02, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + var op3 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.03, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + + // predict the reflectance at each specified optical properties for the given s-d separation + var rOfRho1 = solver.ROfRho(op1, rhos); + var rOfRho2 = solver.ROfRho(op2, rhos); + var rOfRho3 = solver.ROfRho(op3, rhos); + + // Plot log(reflectance) as a function of source-detector separation at each set of optical properties + var chart = Chart.Combine( + new[] { + LineChart(rhos, rOfRho1.Select(r => Math.Log(r)).ToArray(), xLabel: "rho [mm]", yLabel: $"R(rho) [mm-2]", + title: $"log(R(rho) [mm-2])@ mua1={op1[0].RegionOP.Mua:F3}"), + LineChart(rhos, rOfRho2.Select(r => Math.Log(r)).ToArray(), xLabel: "rho [mm]", yLabel: $"R(rho) [mm-2]", + title: $"log(R(rho) [mm-2])@ mua1={op2[0].RegionOP.Mua:F3}"), + LineChart(rhos, rOfRho3.Select(r => Math.Log(r)).ToArray(), xLabel: "rho [mm]", yLabel: $"R(rho) [mm-2]", + title: $"log(R(rho) [mm-2])@ mua1={op3[0].RegionOP.Mua:F3}") + } + ); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo20ROfFxTwoLayerMultiOpProp.cs b/src/Vts.Scripting/ForwardSolvers/Demo20ROfFxTwoLayerMultiOpProp.cs new file mode 100644 index 000000000..3279a255c --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo20ROfFxTwoLayerMultiOpProp.cs @@ -0,0 +1,60 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of spatial frequency +/// for multiple sets of optical properties using a two-layer forward solver +/// +internal class Demo20ROfFxTwoLayerMultiOpProp : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 20: predict R(fx) based on a standard diffusion approximation solution to the time-dependent RTE + // for multiple sets of optical properties using a two-layer forward solver + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new TwoLayerSDAForwardSolver(); + var fxs = new DoubleRange(start: 0, stop: 0.5, number: 51).AsEnumerable().ToArray(); // range of spatial frequencies in 1/mm + var op1 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + var op2 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.02, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + var op3 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.03, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + + // predict the reflectance versus spatial frequency at each specified optical properties + var rOfFx1 = solver.ROfFx(op1, fxs); + var rOfFx2 = solver.ROfFx(op2, fxs); + var rOfFx3 = solver.ROfFx(op3, fxs); + + // Plot reflectance as a function of spatial frequencies at each set of optical properties + var chart = Chart.Combine( + new[] { + LineChart(fxs, rOfFx1, xLabel: "fx [mm-1]", yLabel: $"R(fx) [unitless]", + title: $"R(fx) [unitless] @ mua1={op1[0].RegionOP.Mua:F3}"), + LineChart(fxs, rOfFx2, xLabel: "fx [mm-1]", yLabel: $"R(fx) [unitless]", + title: $"R(fx) [unitless] @ mua1={op2[0].RegionOP.Mua:F3}"), + LineChart(fxs, rOfFx3, xLabel: "fx [mm-1]", yLabel: $"R(fx) [unitless]", + title: $"R(fx) [unitless] @ mua1={op3[0].RegionOP.Mua:F3}") + } + ); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo21ROfRhoAndTimeTwoLayerMultiOpProp.cs b/src/Vts.Scripting/ForwardSolvers/Demo21ROfRhoAndTimeTwoLayerMultiOpProp.cs new file mode 100644 index 000000000..c6aacbe92 --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo21ROfRhoAndTimeTwoLayerMultiOpProp.cs @@ -0,0 +1,87 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of time at a fixed source-detector location +/// for multiple sets of optical properties using a two-layer forward solver +/// +internal class Demo21ROfRhoAndFtTwoLayerMultiOpProp : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 21: predict R(t) at a fixed source-detector location based on a standard diffusion approximation solution + // to the time-dependent RTE for multiple sets of optical properties using a two-layer forward solver + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new TwoLayerSDAForwardSolver(); + var ts = new DoubleRange(start: 0, stop: 0.5, number: 51).AsEnumerable().ToArray(); // range of times in 1/mm + + // create an array of chromophore absorbers, each with a given concentrations + var chromophores = new IChromophoreAbsorber[] + { + new ChromophoreAbsorber(ChromophoreType.HbO2, 70), // molar concentration + new ChromophoreAbsorber(ChromophoreType.Hb, 30), // molar concentration + new ChromophoreAbsorber(ChromophoreType.H2O, 0.8), // fractional concentration + }; + + // construct a scatterer + var scatterer = new PowerLawScatterer(a: 1.2, b: 1.42); + + // compose a tissue using the chromophores and scatterer + var tissue = new Tissue(chromophores, scatterer, "", n: 1.4); + + // predict the bulk tissue's NIR optical properties + var wavelengths = new DoubleRange(start: 650, stop: 850, number: 3).AsEnumerable().ToArray(); // range of wavelengths in nm + var opsBottomLayer = tissue.GetOpticalProperties(wavelengths); + + // perturb the top layer's mua by a multiplicative factor + var muaPerturbationFactor = 1.1; + var opsTopLayer = opsBottomLayer + .Select(op => new OpticalProperties(mua: op.Mua * muaPerturbationFactor, musp: op.Musp, g: op.G, n: op.N)) + .ToArray(); + + // create the multilayer optical properties with the perturbed top layer values + var op1 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: opsTopLayer[0]), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: opsBottomLayer[0]), + }; + var op2 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: opsTopLayer[1]), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: opsBottomLayer[1]), + }; + var op3 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: opsTopLayer[2]), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: opsBottomLayer[2]), + }; + + // predict reflectance versus time at each specified optical properties for the given s-d separation + var rho = 10; // s-d separation in mm + var rOfTime1 = solver.ROfRhoAndTime(op1, rho, ts); + var rOfTime2 = solver.ROfRhoAndTime(op2, rho, ts); + var rOfTime3 = solver.ROfRhoAndTime(op3, rho, ts); + + // Plot reflectance as a function of times at each set of optical properties + var chart = Chart.Combine( + new[] { + LineChart(ts, rOfTime1, xLabel: "t [ns]", yLabel: $"R(t) [mm-2*ns-1] @ rho={rho}mm", + title: $"R(t) [mm-2*s-1] @ rho={rho}mm (mua1={op1[0].RegionOP.Mua:F3}, mua2={op1[1].RegionOP.Mua:F3})"), + LineChart(ts, rOfTime2, xLabel: "t [ns]", yLabel: $"R(t) [mm-2*ns-1] @ rho={rho}mm", + title: $"R(t) [mm-2*s-1] @ rho={rho}mm (mua1={op2[0].RegionOP.Mua:F3}, mua2={op2[1].RegionOP.Mua:F3})"), + LineChart(ts, rOfTime3, xLabel: "t [ns]", yLabel: $"R(t) [mm-2*ns-1] @ rho= {rho}mm", + title: $"R(t) [mm-2*s-1] @ rho={rho}mm (mua1={op3[0].RegionOP.Mua:F3}, mua2={op3[1].RegionOP.Mua:F3})") + } + ); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ForwardSolvers/Demo22ROfRhoAndFtTwoLayerMultiOpProp.cs b/src/Vts.Scripting/ForwardSolvers/Demo22ROfRhoAndFtTwoLayerMultiOpProp.cs new file mode 100644 index 000000000..fd341d27c --- /dev/null +++ b/src/Vts.Scripting/ForwardSolvers/Demo22ROfRhoAndFtTwoLayerMultiOpProp.cs @@ -0,0 +1,94 @@ +namespace Vts.Scripting.ForwardSolvers; + +/// +/// Class using the Vts.dll library to demonstrate predicting reflectance as a function of time at a fixed source-detector location +/// for multiple sets of optical properties using a two-layer forward solver +/// +internal class Demo22ROfRhoAndFtTwoLayerMultiOpProp : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 22: predict R(ft) at a fixed source-detector location based on a standard diffusion approximation solution + // to the time-dependent RTE for multiple sets of optical properties using a two-layer forward solver + + // Solver type options: + // PointSourceSDA,DistributedGaussianSourceSDA, DistributedPointSourceSDA, + // MonteCarlo(basic scaled), Nurbs(scaled with smoothing and adaptive binning) + var solver = new TwoLayerSDAForwardSolver(); + var fts = new DoubleRange(start: 0, stop: 0.5, number: 51).AsEnumerable().ToArray(); // range of time frequencies in GHz + var op1 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + var op2 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.02, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + var op3 = new IOpticalPropertyRegion[] + { + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.03, musp: 1, g: 0.8, n: 1.4)), + new LayerOpticalPropertyRegion(zRange: new DoubleRange(0, 2, 2), regionOP: new OpticalProperties(mua: 0.01, musp: 1, g: 0.8, n: 1.4)), + }; + + // predict the reflectance at each specified optical properties for the given s-d separation + var rho = 10; // s-d separation in mm + var rOfFt1 = solver.ROfRhoAndFt(op1, rho, fts); + var rOfFt2 = solver.ROfRhoAndFt(op2, rho, fts); + var rOfFt3 = solver.ROfRhoAndFt(op3, rho, fts); + + // Plot reflectance amplitude and phase as a function of times at each set of optical properties + var magnitudeScale = 1E4; + var amplitudeChart = Chart.Combine( + new[] { + LineChart(fts, rOfFt1.Select(r => r.Magnitude * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"R(ft) [mm-2] * {magnitudeScale:E0}", + title: $"|R(ft)| [mm-2] @ rho={rho}mm (mua1={op1[0].RegionOP.Mua:F3}, mua2={op1[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt2.Select(r => r.Magnitude * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"R(ft) [mm-2] * {magnitudeScale:E0}", + title: $"|R(ft)| [mm-2] @ rho={rho}mm (mua1={op2[0].RegionOP.Mua:F3}, mua2={op2[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt3.Select(r => r.Magnitude * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"R(ft) [mm-2] * {magnitudeScale:E0}", + title: $"|R(ft)| [mm-2] @ rho={rho}mm (mua1={op3[0].RegionOP.Mua:F3}, mua2={op3[1].RegionOP.Mua:F3})") + } + ); + var phaseChart = Chart.Combine( + new[] { + LineChart(fts, rOfFt1.Select(r => -r.Phase).ToArray(), xLabel: "ft [GHz]", yLabel: $"Φ(R(ft)) [rad] @ rho={rho}mm", + title: $"Φ(R(ft)) [rad] @ rho={rho}mm (mua1={op1[0].RegionOP.Mua:F3}, mua2={op1[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt2.Select(r => -r.Phase).ToArray(), xLabel: "ft [GHz]", yLabel: $"Φ(R(ft)) [rad] @ rho={rho}mm", + title: $"Φ(R(ft)) [rad] @ rho={rho}mm (mua1={op2[0].RegionOP.Mua:F3}, mua2={op2[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt3.Select(r => -r.Phase).ToArray(), xLabel: "ft [GHz]", yLabel: $"Φ(R(ft)) [rad] @ rho={rho}mm", + title: $"Φ(R(ft)) [rad] @ rho={rho}mm (mua1={op3[0].RegionOP.Mua:F3}, mua2={op3[1].RegionOP.Mua:F3})") + } + ); + var realChart = Chart.Combine( + new[] { + LineChart(fts, rOfFt1.Select(r => r.Real * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"Re(R(ft)) [mm-2]*{magnitudeScale:E0}", + title: $"Re(R(ft)) [mm-2] @ rho={rho}mm (mua1={op1[0].RegionOP.Mua:F3}, mua2={op1[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt2.Select(r => r.Real * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"Re(R(ft)) [mm-2]*{magnitudeScale:E0}", + title: $"Re(R(ft)) [mm-2] @ rho={rho}mm (mua1={op2[0].RegionOP.Mua:F3}, mua2={op2[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt3.Select(r => r.Real * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"Re(R(ft)) [mm-2]*{magnitudeScale:E0}", + title: $"Re(R(ft)) [mm-2] @ rho={rho}mm (mua1={op3[0].RegionOP.Mua:F3}, mua2={op3[1].RegionOP.Mua:F3})") + } + ); + var imagChart = Chart.Combine( + new[] { + LineChart(fts, rOfFt1.Select(r => r.Imaginary * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"Imag(R(ft)) [mm-2]*{magnitudeScale:E0}", + title: $"Imag(R(ft)) [mm-2] @ rho={rho}mm (mua1={op1[0].RegionOP.Mua:F3}, mua2={op1[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt2.Select(r => r.Imaginary * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"Imag(R(ft)) [mm-2]*{magnitudeScale:E0}", + title: $"Imag(R(ft)) [mm-2] @ rho={rho}mm (mua1={op2[0].RegionOP.Mua:F3}, mua2={op2[1].RegionOP.Mua:F3})"), + LineChart(fts, rOfFt3.Select(r => r.Imaginary * magnitudeScale).ToArray(), xLabel: "ft [GHz]", yLabel: $"Imag(R(ft)) [mm-2]*{magnitudeScale:E0}", + title: $"Imag(R(ft)) [mm-2] @ rho={rho}mm (mua1={op3[0].RegionOP.Mua:F3}, mua2={op3[1].RegionOP.Mua:F3})") + } + ); + var realImagChart = Chart.Grid(new[] { realChart, imagChart }, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + var ampPhaseChart = Chart.Grid(new[] { amplitudeChart, phaseChart }, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + if (showPlots) + { + realImagChart.Show(); + ampPhaseChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/GlobalUsings.cs b/src/Vts.Scripting/GlobalUsings.cs new file mode 100644 index 000000000..514d943b0 --- /dev/null +++ b/src/Vts.Scripting/GlobalUsings.cs @@ -0,0 +1,17 @@ +global using Vts.Common; +global using Vts.Extensions; +global using Vts.Modeling.Optimizers; +global using Vts.Modeling.ForwardSolvers; +global using Vts.SpectralMapping; +global using Vts.Factories; +global using Vts.MonteCarlo; +global using Vts.MonteCarlo.Sources; +global using Vts.MonteCarlo.Tissues; +global using Vts.MonteCarlo.Detectors; +global using Vts.MonteCarlo.Factories; +global using Vts.MonteCarlo.PhotonData; +global using Vts.MonteCarlo.PostProcessing; +global using static Vts.Scripting.ScriptHelper; + +global using System.Numerics; +global using Plotly.NET.CSharp; \ No newline at end of file diff --git a/src/Vts.Scripting/IDemoScript.cs b/src/Vts.Scripting/IDemoScript.cs new file mode 100644 index 000000000..65d2f116e --- /dev/null +++ b/src/Vts.Scripting/IDemoScript.cs @@ -0,0 +1,12 @@ +namespace Vts.Scripting; + +/// +/// Interface to constrain all demo scripts to a uniform signature +/// +interface IDemoScript +{ + /// + /// The one required static method for all demo scripts + /// + // static abstract void RunDemo(bool showPlots = true); // todo: fix this in a future branch +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo01ROfRhoSimple.cs b/src/Vts.Scripting/MonteCarlo/Demo01ROfRhoSimple.cs new file mode 100644 index 000000000..0c4ecc5a2 --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo01ROfRhoSimple.cs @@ -0,0 +1,46 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate performing a simple Monte Carlo simulation +/// +internal class Demo01ROfRhoSimple : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 01: run a simple Monte Carlo simulation with 1000 photons. + // Notes: + // - default source is a point source beam normally incident at the origin + // - default tissue is a 100mm thick slab with air-tissue boundary and optical properties: mua: 0.01, musp: 1.0, g: 0.8, n:1.4 + + // create a SimulationInput object to define the simulation + var detectorRange = new DoubleRange(start: 0, stop: 40, number: 201); + var simulationInput = new SimulationInput + { + // specify the number of photons to run + N = 1000, + + // define a single R(rho) detector by the endpoints of rho bins + DetectorInputs = new List { new ROfRhoDetectorInput { Rho = detectorRange, Name = "ROfRho" } }, // name can be whatever you want + }; + + // create the simulation + var simulation = new MonteCarloSimulation(simulationInput); + + // run the simulation + var simulationOutput = simulation.Run(); + + // plot the results using Plotly.NET + var detectorResults = (ROfRhoDetector)simulationOutput.ResultsDictionary["ROfRho"]; + var logReflectance = detectorResults.Mean.Select(r => Math.Log(r)).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (detectorRange.GetMidpoints(), "ρ [mm]", "log(R(ρ)) [mm-2]"); + var chart = LineChart(detectorMidpoints, logReflectance, xLabel, yLabel, title: "log(R(ρ)) [mm-2]"); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo02DAWvsCAW.cs b/src/Vts.Scripting/MonteCarlo/Demo02DAWvsCAW.cs new file mode 100644 index 000000000..7eee532af --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo02DAWvsCAW.cs @@ -0,0 +1,72 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate comparing Monte Carlo simulations with different absorption weighting types +/// +internal class Demo02DawVsCaw : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 02: run Monte Carlo simulations for two absorption weighting types + // with 1000 photons each and compare computation time + // Notes: + // - default source is a point source beam normally incident at the origin + // - default tissue is a 100mm thick slab with air-tissue boundary and optical properties: mua: 0.01, musp: 1.0, g: 0.8, n:1.4 + + // define some shared values between the two simulations + var numPhotons = 1000; + var detectorRange = new DoubleRange(start: 0, stop: 40, number: 201); + var detectorInput = new ROfRhoDetectorInput { Rho = detectorRange, Name = "ROfRho" }; // name can be whatever you want + + // create SimulationInput objects to define the two simulations + var simulationInput1 = new SimulationInput + { + N = numPhotons, + DetectorInputs = new List { detectorInput }, + OutputName = "MonteCarlo02ROfRho-DAW", + Options = new SimulationOptions + { + AbsorptionWeightingType = AbsorptionWeightingType.Discrete, // discrete absorption weighting + }, + }; + + var simulationInput2 = new SimulationInput + { + N = numPhotons, + DetectorInputs = new List { detectorInput }, + OutputName = "MonteCarlo02ROfRho-CAW", + Options = new SimulationOptions + { + AbsorptionWeightingType = AbsorptionWeightingType.Continuous, // continuous absorption weighting + } + }; + + // create the simulations + var simulation1 = new MonteCarloSimulation(simulationInput1); + var simulation2 = new MonteCarloSimulation(simulationInput2); + + // run the simulations + var simulation1Output = simulation1.Run(); + var simulation2Output = simulation2.Run(); + + // plot and compare the results using Plotly.NET + var detectorResults1 = (ROfRhoDetector)simulation1Output.ResultsDictionary[detectorInput.Name]; + var detectorResults2 = (ROfRhoDetector)simulation2Output.ResultsDictionary[detectorInput.Name]; + var logReflectance1 = detectorResults1.Mean.Select(r => Math.Log(r)).ToArray(); + var logReflectance2 = detectorResults2.Mean.Select(r => Math.Log(r)).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (detectorRange.GetMidpoints(), "ρ [mm]", "log(R(ρ)) [mm-2]"); + var chart = Chart.Combine(new[] + { + LineChart(detectorMidpoints, logReflectance1, xLabel, yLabel, title: "log(R(ρ)) [mm-2] - CAW"), + LineChart(detectorMidpoints, logReflectance2, xLabel, yLabel, title: "log(R(ρ)) [mm-2] - DAW") + }); // show both charts together + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo03ROfRhoFullCustomization.cs b/src/Vts.Scripting/MonteCarlo/Demo03ROfRhoFullCustomization.cs new file mode 100644 index 000000000..46d6ed89e --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo03ROfRhoFullCustomization.cs @@ -0,0 +1,82 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate performing a Monte Carlo simulation with full customization +/// +internal class Demo03ROfRhoFullCustomization : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 03: run a Monte Carlo simulation with a fully-customized input + // (values used here are the class defaults) + + var detectorRange = new DoubleRange(start: 0, stop: 40, number: 201); + + // create a SimulationInput object to define the simulation + var simulationInput = new SimulationInput + { + // specify the number of photons + N = 1000, + + // define a point source beam normally incident at the origin + SourceInput = new DirectionalPointSourceInput + { + SourceType = "DirectionalPoint", + PointLocation = new(x: 0, y: 0, z: 0), + Direction = new(ux: 0, uy: 0, uz: 1), + InitialTissueRegionIndex = 0 + }, + + // define a semi-infinite slab tissue geometry with air-tissue boundary (a bottom air layer is necessary) + TissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 100), // tissue "z" range ("semi-infinite" slab, 100mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(100, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }, + + // define a single R(rho) detector by the endpoints of rho bins + DetectorInputs = new IDetectorInput[] + { + new ROfRhoDetectorInput { Rho = detectorRange, TallySecondMoment = true, Name = "ROfRho" } // name can be whatever you want + }, + + // specify all simulation options + Options = new SimulationOptions + { + Seed = 0, // -1 will generate a random seed + AbsorptionWeightingType = AbsorptionWeightingType.Discrete, + PhaseFunctionType = PhaseFunctionType.HenyeyGreenstein + } + }; + + // create the simulation + var simulation = new MonteCarloSimulation(simulationInput); + + // run the simulation + var simulationOutput = simulation.Run(); + + // plot the results using Plotly.NET + var detectorResults = (ROfRhoDetector)simulationOutput.ResultsDictionary["ROfRho"]; + var logReflectance = detectorResults.Mean.Select(r => Math.Log(r)).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (detectorRange.GetMidpoints(), "ρ [mm]", "log(R(ρ)) [mm-2]"); + var chart = LineChart(detectorMidpoints, logReflectance, xLabel, yLabel, title: "log(R(ρ)) [mm-2]"); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo04N1000vsN100.cs b/src/Vts.Scripting/MonteCarlo/Demo04N1000vsN100.cs new file mode 100644 index 000000000..3d4935a8a --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo04N1000vsN100.cs @@ -0,0 +1,59 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate comparing two Monte Carlo simulations with different photon counts +/// +internal class Demo04N1000VsN100 : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 04: run a list of Monte Carlo simulations + // create a list of two default SimulationInput with different numbers of photons + // Notes: + // - default source is a point source beam normally incident at the origin + // - default tissue is a 100mm thick slab with air-tissue boundary + + // define some shared values between the two simulations + var detectorRange = new DoubleRange(start: 0, stop: 40, number: 201); + var detectorInput = new ROfRhoDetectorInput { Rho = detectorRange, Name = "ROfRho" }; + + // create two simulations, one with n=1000 and one with n=100 + // note, we are creating a SimulationInput object here in-line + var simulation1 = new MonteCarloSimulation(input: new SimulationInput + { + N = 1000, + DetectorInputs = new IDetectorInput[] { detectorInput }, + OutputName = "MonteCarlo04ROfRho-n1000" + }); + var simulation2 = new MonteCarloSimulation(input: new SimulationInput + { + N = 100, + DetectorInputs = new IDetectorInput[] { detectorInput }, + OutputName = "MonteCarlo04ROfRho-n100" + }); + + // run the simulations + var simulation1Output = simulation1.Run(); + var simulation2Output = simulation2.Run(); + + // plot and compare the results using Plotly.NET + var detectorResults1 = (ROfRhoDetector)simulation1Output.ResultsDictionary[detectorInput.Name]; + var detectorResults2 = (ROfRhoDetector)simulation2Output.ResultsDictionary[detectorInput.Name]; + var logReflectance1 = detectorResults1.Mean.Select(r => Math.Log(r)).ToArray(); + var logReflectance2 = detectorResults2.Mean.Select(r => Math.Log(r)).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (detectorRange.GetMidpoints(), "ρ [mm]", "log(R(ρ)) [mm-2]"); + var chart = Chart.Combine(new[] + { + LineChart(detectorMidpoints, logReflectance1, xLabel, yLabel, title: "log(R(ρ)) [mm-2] - n=1000"), + LineChart(detectorMidpoints, logReflectance2, xLabel, yLabel, title: "log(R(ρ)) [mm-2] - n=100") + }); // show both charts together + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo05PostProcessor.cs b/src/Vts.Scripting/MonteCarlo/Demo05PostProcessor.cs new file mode 100644 index 000000000..9ca05d9d3 --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo05PostProcessor.cs @@ -0,0 +1,74 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate using the Monte Carlo post-processor +/// +internal class Demo05PostProcessor : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 05: run a Monte Carlo simulation with post-processing enabled. First run a simulation, + // then post-process the generated database and compare on-the-fly results with post-processed results + // Notes: + // - the on-the-fly ROfRhoDetectorInput is not necessary for post-processing, but is included here for comparison + // - default source is a point source beam normally incident at the origin + // - default tissue is a 100mm thick slab with air-tissue boundary + + // create a SimulationInput object to define the simulation + var detectorRange = new DoubleRange(start: 0, stop: 40, number: 201); + var detectorInput = new ROfRhoDetectorInput { Rho = detectorRange, Name = "ROfRho" }; + var simulationInput = new SimulationInput + { + N = 1000, + DetectorInputs = new IDetectorInput[] { detectorInput }, + OutputName = "results", + Options = new SimulationOptions + { + Databases = new [] { DatabaseType.DiffuseReflectance } + } + }; + + // create and run the simulation + var simulation = new MonteCarloSimulation(input: simulationInput); + var simulationOutput = simulation.Run(); + + // specify post-processing of generated database + var postProcessorInput = new PostProcessorInput + { + InputFolder = simulationInput.OutputName, + DetectorInputs = new IDetectorInput[] { detectorInput }, + DatabaseSimulationInputFilename = "infile" + }; + + // create and run the post-processor + var photonDatabase = PhotonDatabaseFactory.GetPhotonDatabase( + virtualBoundaryType: VirtualBoundaryType.DiffuseReflectance, + filePath: simulationInput.OutputName); + var postProcessor = new PhotonDatabasePostProcessor( + virtualBoundaryType: VirtualBoundaryType.DiffuseReflectance, + detectorInputs: postProcessorInput.DetectorInputs, + photonDatabase: photonDatabase, + databaseInput: simulationInput); + var postProcessorOutput = postProcessor.Run(); + + // plot and compare the results using Plotly.NET + var onTheFlyDetectorResults = (ROfRhoDetector)simulationOutput.ResultsDictionary[detectorInput.Name]; + var postProcessorDetectorResults = (ROfRhoDetector)postProcessorOutput.ResultsDictionary[detectorInput.Name]; + var logReflectance1 = onTheFlyDetectorResults.Mean.Select(r => Math.Log(r)).ToArray(); + var logReflectance2 = postProcessorDetectorResults.Mean.Select(r => Math.Log(r)).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (detectorRange.GetMidpoints(), "ρ [mm]", "log(R(ρ)) [mm-2]"); + var chart = Chart.Combine(new[] + { + LineChart(detectorMidpoints, logReflectance1, xLabel, yLabel, title: "log(R(ρ)) [mm-2] - on the fly"), + LineChart(detectorMidpoints, logReflectance2, xLabel, yLabel, title: "log(R(ρ)) [mm-2] - post-processor") + }); // show both charts together + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo06pMCPostProcessor.cs b/src/Vts.Scripting/MonteCarlo/Demo06pMCPostProcessor.cs new file mode 100644 index 000000000..de0b076f7 --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo06pMCPostProcessor.cs @@ -0,0 +1,121 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate using the Perturbation Monte Carlo post-processor +/// +internal class Demo06PmcPostProcessor : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 06: run a Monte Carlo simulation with Perturbation Monte Carlo (pMC) post-processing enabled. + // First run a simulation, then post-process the generated database with pMC, varying optical properties post-simulation + // Notes: + // - default source is a point source beam normally incident at the origin + + // create a SimulationInput object to define the simulation + var tissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 100), // tissue "z" range ("semi-infinite" slab, 100mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(100, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }; + var detectorRange = new DoubleRange(start: 0, stop: 40, number: 201); + var simulationInput = new SimulationInput + { + N = 1000, + TissueInput = tissueInput, + DetectorInputs = new IDetectorInput[] { }, // leaving this empty - no on-the-fly detectors needed! + OutputName = "results", + Options = new SimulationOptions + { + Seed = 0, + Databases = new [] { DatabaseType.pMCDiffuseReflectance } + } + }; + + // create and run the simulation, creating the database + var simulation = new MonteCarloSimulation(input: simulationInput); + var simulationOutput = simulation.Run(); + + // specify post-processing of three post-processor detectors with 1x, 0.5x, and 2x the baseline mua (0.01) + var pmcDetectorInput1xMua = new pMCROfRhoDetectorInput + { + Rho = detectorRange, + Name = "ROfRho", + PerturbedOps = new OpticalProperties[] + { + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0), // air optical properties + new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4), // tissue optical properties (baseline) + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0) // air optical properties + }, + PerturbedRegionsIndices = new[] { 1 } // only perturbing the tissue OPs + }; + var pmcDetectorInput0p5xMua = new pMCROfRhoDetectorInput + { + Rho = detectorRange, + Name = "ROfRhoMuaHalf", + PerturbedOps = new OpticalProperties[] + { + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0), // air optical properties + new(mua: 0.005, musp: 1.0, g: 0.8, n: 1.4), // tissue optical properties (half baseline) + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0) // air optical properties + }, + PerturbedRegionsIndices = new[] { 1 } // only perturbing the tissue OPs + }; + var pmcDetectorInput2xMua = new pMCROfRhoDetectorInput + { + Rho = detectorRange, + Name = "ROfRhoMuaDouble", + PerturbedOps = new OpticalProperties[] + { + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0), // air optical properties + new(mua: 0.02, musp: 1.0, g: 0.8, n: 1.4), // tissue optical properties (double baseline) + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0) // air optical properties + }, + PerturbedRegionsIndices = new[] { 1 } // only perturbing the tissue OPs + }; + + // create and run the post-processor + var photonDatabase = PhotonDatabaseFactory.GetpMCDatabase( + virtualBoundaryType: VirtualBoundaryType.pMCDiffuseReflectance, + filePath: simulationOutput.Input.OutputName); + var postProcessor = new PhotonDatabasePostProcessor( + virtualBoundaryType: VirtualBoundaryType.pMCDiffuseReflectance, + detectorInputs: new IDetectorInput[] { pmcDetectorInput1xMua, pmcDetectorInput0p5xMua, pmcDetectorInput2xMua }, + database: photonDatabase, + databaseInput: simulationOutput.Input); + var postProcessorOutput = postProcessor.Run(); + + // plot and compare the results using Plotly.NET + var postProcessorDetectorResults1xMua = (pMCROfRhoDetector)postProcessorOutput.ResultsDictionary[pmcDetectorInput1xMua.Name]; + var postProcessorDetectorResults0p5xMua = (pMCROfRhoDetector)postProcessorOutput.ResultsDictionary[pmcDetectorInput0p5xMua.Name]; + var postProcessorDetectorResults2xMua = (pMCROfRhoDetector)postProcessorOutput.ResultsDictionary[pmcDetectorInput2xMua.Name]; + var logReflectance1 = postProcessorDetectorResults1xMua.Mean.Select(r => Math.Log(r)).ToArray(); + var logReflectance2 = postProcessorDetectorResults0p5xMua.Mean.Select(r => Math.Log(r)).ToArray(); + var logReflectance3 = postProcessorDetectorResults2xMua.Mean.Select(r => Math.Log(r)).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (detectorRange.GetMidpoints(), "ρ [mm]", "log(R(ρ)) [mm-2]"); + var chart = Chart.Combine(new[] + { + LineChart(detectorMidpoints, logReflectance1, xLabel, yLabel, title: $"log(R(ρ)) [mm-2] - 1.0x baseline (mua={0.01:F3}/mm)"), + LineChart(detectorMidpoints, logReflectance2, xLabel, yLabel, title: $"log(R(ρ)) [mm-2] - 0.5x baseline (mua={0.005:F3}/mm)"), + LineChart(detectorMidpoints, logReflectance3, xLabel, yLabel, title: $"log(R(ρ)) [mm-2] - 2.0x baseline (mua={0.02:F3}/mm)") + }); // show all three charts together + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo07pMCInversion.cs b/src/Vts.Scripting/MonteCarlo/Demo07pMCInversion.cs new file mode 100644 index 000000000..edfd1717f --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo07pMCInversion.cs @@ -0,0 +1,664 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate using the Perturbation Monte Carlo post-processor to calculate optical properties (i.e. "inversion") +/// +internal class Demo07PmcInversion : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 07: run a Monte Carlo simulation with Perturbation Monte Carlo post-processor to calculate + // optical properties (i.e. "inversion") based on fitting to measured data generated using Nurbs + // Notes: + // - default source is a point source beam normally incident at the origin + // - convergence to measured data optical properties affected by: + // - number of photons launched in baseline simulation, N + // - placement and number of rho + // - distance of initial guess from actual + // - normalization of chi2 + // - weighting of chi2 via the standard deviation provided + // - optimization options selected + + // create a SimulationInput object to define the simulation + var tissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 100), // tissue "z" range ("semi-infinite" slab, 100mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(100, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }; + var simulationInput = new SimulationInput + { + N = 1000, + TissueInput = tissueInput, + DetectorInputs = new IDetectorInput[] { }, // leaving this empty - no on-the-fly detectors needed! + OutputName = "results", + Options = new SimulationOptions + { + //Seed = 0, + AbsorptionWeightingType = AbsorptionWeightingType.Discrete, + Databases = new [] { DatabaseType.pMCDiffuseReflectance } + } + }; + + // create and run the simulation, generating the pMC database + var simulation = new MonteCarloSimulation(input: simulationInput); + var simulationOutput = simulation.Run(); + + // Create some measurements, based on a Nurbs-based White Monte Carlo forward solver + var measuredOPs = new OpticalProperties(mua: 0.04, musp: 0.95, g: 0.8, n: 1.4); + var detectorRange = new DoubleRange(start: 0, stop: 6, number: 7); + var detectorMidpoints = detectorRange.GetMidpoints(); + var measurementForwardSolver = new NurbsForwardSolver(); + var measuredData = measurementForwardSolver.ROfRho(measuredOPs, detectorMidpoints); + + // Specify initial guess for optimization equal to the original pMC simulation baseline values + var baselineOps = tissueInput.Regions[1].RegionOP; + + // Create an ad-hoc forward solver based on pMC prediction (see implementation below; note: implemented for ROfRho only) + var pmcForwardSolver = new PmcForwardSolver(detectorRange, simulationOutput.Input); + + // Run the inversion, based on Levenberg-Marquardt optimization. + // Note: chi-squared weighting in the inversion is based on standard deviation of measured data, except for + // the last point, which is set to infinity to fully-deweight it. This is because the last point contains + // all photon counts for that bin and beyond, and we therefore don't want to include it in the fit + var initialGuessOPsAndRhoAxis = new object[] { new[] { baselineOps }, detectorMidpoints }; // "extra" (constant) data for the forward model calls + var solution = ComputationFactory.SolveInverse( + forwardSolver: pmcForwardSolver, + optimizer: new MPFitLevenbergMarquardtOptimizer(), + solutionDomainType: SolutionDomainType.ROfRho, + dependentValues: measuredData, + standardDeviationValues: measuredData[..^1].Append(double.PositiveInfinity).ToArray(), + inverseFitType: InverseFitType.MuaMusp, + initialGuessOPsAndRhoAxis + ); + // convert the returned double[] {mua, musp} array to an OpticalProperties object + var fitOps = new OpticalProperties(solution[0], solution[1], baselineOps.G, baselineOps.N); + + // For illustration purposes, also run the pMC solver + var postProcessorDetectorResultsInitialGuess = pmcForwardSolver.ROfRho(baselineOps, detectorMidpoints); + var postProcessorDetectorResultsFitValues = pmcForwardSolver.ROfRho(fitOps, detectorMidpoints); + + // plot and compare the results using Plotly.NET + var logReflectance1 = postProcessorDetectorResultsInitialGuess.Select(r => Math.Log(r)).ToArray(); + var logReflectance2 = postProcessorDetectorResultsFitValues.Select(r => Math.Log(r)).ToArray(); + var (xLabel, yLabel) = ("ρ [mm]", "log(R(ρ)) [mm-2]"); + var chart = Chart.Combine(new[] + { + ScatterChart(detectorMidpoints, logReflectance2, xLabel, yLabel, title: "Measured Data"), + LineChart(detectorMidpoints, logReflectance1, xLabel, yLabel, title: $"Initial guess (mua={baselineOps.Mua:F3}/mm, musp={baselineOps.Musp:F3}/mm)"), + LineChart(detectorMidpoints, logReflectance2, xLabel, yLabel, title: $"Converged result (mua={fitOps.Mua:F3}/mm, musp={fitOps.Musp:F3}/mm)") + }); // show all three charts together + + if (showPlots) + { + chart.Show(); + } + } + + /// + /// The PmcForwardSolver class is a private class that implements the IForwardSolver interface. + /// It is used for performing forward simulations using a perturbation Monte Carlo (pMC) approach. + /// The class uses a specified detector range, simulation input parameters, and a pMC database + /// to calculate reflectance as a function of radial distance for given sets of optical properties. + /// + private class PmcForwardSolver : IForwardSolver + { + private readonly DoubleRange _detectorRange; + private readonly SimulationInput _simulationInput; + private readonly pMCDatabase _photonDatabase; + + /// + /// Initializes a new instance of the PmcForwardSolver class using the specified detector range, simulation input, and photon database. + /// + /// The range of the detector for the simulation. + /// The input parameters for the simulation. + /// The pMC (probability Monte Carlo) database used for the simulation. + /// + /// This is a private constructor used internally by the class for initialization. + /// + private PmcForwardSolver(DoubleRange detectorRange, SimulationInput simulationInput, pMCDatabase photonDatabase) + { + _detectorRange = detectorRange; + _simulationInput = simulationInput; + _photonDatabase = photonDatabase; + } + + /// + /// Initializes a new instance of the PmcForwardSolver class using the specified detector range and simulation input. + /// + /// The range of the detector for the simulation. + /// The input parameters for the simulation. + /// + /// This constructor uses the output name from the simulation input to create a pMC (perturbation Monte Carlo) database. + /// The pMC database is retrieved using the PhotonDatabaseFactory's GetpMCDatabase method with the virtual boundary type set to pMCDiffuseReflectance. + /// + public PmcForwardSolver(DoubleRange detectorRange, SimulationInput simulationInput) + : this(detectorRange, simulationInput, PhotonDatabaseFactory.GetpMCDatabase( + virtualBoundaryType: VirtualBoundaryType.pMCDiffuseReflectance, + filePath: simulationInput.OutputName)) + { + } + + /// + /// Calculates the reflectance as a function of radial distance (rho) for a given set of optical properties. + /// + /// The optical properties of the medium. + /// An array of radial distances at which the reflectance is to be calculated. + /// An array of reflectance values corresponding to each radial distance in the input array. + public double[] ROfRho(OpticalProperties op, double[] rhos) + { + // create and run the post-processor + var pmcDetectorInput = new pMCROfRhoDetectorInput + { + Rho = _detectorRange, + Name = "ROfRho-OnTheFly", + PerturbedOps = new OpticalProperties[] + { + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0), // air optical properties + new(mua: op.Mua, musp: op.Musp, g: op.G, n: op.N), // tissue optical properties (perturbed) + new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0) // air optical properties + }, + PerturbedRegionsIndices = new[] { 1 } // only perturbing the tissue OPs + }; + var postProcessor = new PhotonDatabasePostProcessor( + virtualBoundaryType: VirtualBoundaryType.pMCDiffuseReflectance, + detectorInputs: new IDetectorInput[] { pmcDetectorInput }, + database: _photonDatabase, + databaseInput: _simulationInput); + var postProcessorOutput = postProcessor.Run(); + var postProcessorDetectorResults = (pMCROfRhoDetector)postProcessorOutput.ResultsDictionary[pmcDetectorInput.Name]; + return postProcessorDetectorResults.Mean; + } + + /// + /// Calculates the reflectance as a function of radial distance (rho) for each set of optical properties provided. + /// + /// An array of optical properties of the medium. + /// An array of radial distances at which the reflectance is to be calculated. + /// An array of reflectance values corresponding to each set of optical properties and each radial distance in the input arrays. + public double[] ROfRho(OpticalProperties[] ops, double[] rhos) + { + return ops.SelectMany(op => ROfRho(op, rhos)).ToArray(); + } + + #region not implemented + public double BeamDiameter { get => throw new NotImplementedException(); set => throw new NotImplementedException(); } + + public event System.ComponentModel.PropertyChangedEventHandler? PropertyChanged; + + public double ROfRho(IOpticalPropertyRegion[] regions, double rho) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfRho(IEnumerable regions, IEnumerable rhos) + { + throw new NotImplementedException(); + } + + public double[] ROfRho(IOpticalPropertyRegion[] regions, double[] rhos) + { + throw new NotImplementedException(); + } + + public double[] ROfRho(IOpticalPropertyRegion[][] regions, double[] rhos) + { + throw new NotImplementedException(); + } + + public double ROfRhoAndTime(IOpticalPropertyRegion[] regions, double rho, double time) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfRhoAndTime(IEnumerable regions, IEnumerable rhos, IEnumerable times) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(IOpticalPropertyRegion[] regions, double[] rhos, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(IOpticalPropertyRegion[] regions, double[] rhos, double time) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(IOpticalPropertyRegion[] regions, double rho, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(IOpticalPropertyRegion[][] regions, double[] rhos, double[] times) + { + throw new NotImplementedException(); + } + + public Complex ROfRhoAndFt(IOpticalPropertyRegion[] regions, double rho, double ft) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfRhoAndFt(IEnumerable regions, IEnumerable rhos, IEnumerable fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(IOpticalPropertyRegion[] regions, double rho, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(IOpticalPropertyRegion[] regions, double[] rhos, double ft) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(IOpticalPropertyRegion[] regions, double[] rhos, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(IOpticalPropertyRegion[][] regions, double[] rhos, double[] fts) + { + throw new NotImplementedException(); + } + + public double ROfFx(IOpticalPropertyRegion[] regions, double fx) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfFx(IEnumerable regions, IEnumerable fxs) + { + throw new NotImplementedException(); + } + + public double[] ROfFx(IOpticalPropertyRegion[] regions, double[] fxs) + { + throw new NotImplementedException(); + } + + public double[] ROfFx(IOpticalPropertyRegion[][] regions, double[] fxs) + { + throw new NotImplementedException(); + } + + public double ROfFxAndTime(IOpticalPropertyRegion[] regions, double fx, double time) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfFxAndTime(IEnumerable regions, IEnumerable fxs, IEnumerable times) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(IOpticalPropertyRegion[] regions, double[] fxs, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(IOpticalPropertyRegion[] regions, double fx, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(IOpticalPropertyRegion[] regions, double[] fxs, double time) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(IOpticalPropertyRegion[][] regions, double[] fxs, double[] times) + { + throw new NotImplementedException(); + } + + public Complex ROfFxAndFt(IOpticalPropertyRegion[] regions, double fx, double ft) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfFxAndFt(IEnumerable regions, IEnumerable fxs, IEnumerable fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(IOpticalPropertyRegion[] regions, double fx, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(IOpticalPropertyRegion[] regions, double[] fxs, double ft) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(IOpticalPropertyRegion[] regions, double[] fxs, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(IOpticalPropertyRegion[][] regions, double[] fxs, double[] fts) + { + throw new NotImplementedException(); + } + + public double[] FluenceOfRhoAndZ(IOpticalPropertyRegion[][] regions, double[] rhos, double[] zs) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfRhoAndZ(IEnumerable regions, IEnumerable rhos, IEnumerable zs) + { + throw new NotImplementedException(); + } + + public Complex[] FluenceOfRhoAndZAndFt(IOpticalPropertyRegion[][] regions, double[] rhos, double[] zs, double[] fts) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfRhoAndZAndFt(IEnumerable regions, IEnumerable rhos, IEnumerable zs, IEnumerable fts) + { + throw new NotImplementedException(); + } + + public double ROfRho(OpticalProperties op, double rho) + { + throw new NotImplementedException(); + } + + public double ROfRhoAndTime(OpticalProperties op, double rho, double time) + { + throw new NotImplementedException(); + } + + public Complex ROfRhoAndFt(OpticalProperties op, double rho, double ft) + { + throw new NotImplementedException(); + } + + public double ROfFx(OpticalProperties op, double fx) + { + throw new NotImplementedException(); + } + + public double ROfFxAndTime(OpticalProperties op, double fx, double time) + { + throw new NotImplementedException(); + } + + public Complex ROfFxAndFt(OpticalProperties op, double fx, double ft) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfRhoAndZ(IEnumerable ops, IEnumerable rhos, IEnumerable zs) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfRhoAndZAndTime(IEnumerable ops, IEnumerable rhos, IEnumerable zs, IEnumerable times) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfRhoAndZAndFt(IEnumerable ops, IEnumerable rhos, IEnumerable zs, IEnumerable fts) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfFxAndZ(IEnumerable ops, IEnumerable fxs, IEnumerable zs) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfFxAndZAndTime(IEnumerable ops, IEnumerable fxs, IEnumerable zs, IEnumerable times) + { + throw new NotImplementedException(); + } + + public IEnumerable FluenceOfFxAndZAndFt(IEnumerable ops, IEnumerable fxs, IEnumerable zs, IEnumerable fts) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfRho(IEnumerable ops, IEnumerable rhos) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfRhoAndTime(IEnumerable ops, IEnumerable rhos, IEnumerable times) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfRhoAndFt(IEnumerable ops, IEnumerable rhos, IEnumerable fts) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfFx(IEnumerable ops, IEnumerable fxs) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfFxAndTime(IEnumerable ops, IEnumerable fxs, IEnumerable times) + { + throw new NotImplementedException(); + } + + public IEnumerable ROfFxAndFt(IEnumerable ops, IEnumerable fxs, IEnumerable fts) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(OpticalProperties op, double[] rhos, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(OpticalProperties[] ops, double[] rhos, double[] times) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(OpticalProperties op, double[] rhos, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(OpticalProperties[] ops, double[] rhos, double[] fts) + { + throw new NotImplementedException(); + } + + public double[] ROfFx(OpticalProperties op, double[] fxs) + { + throw new NotImplementedException(); + } + + public double[] ROfFx(OpticalProperties[] ops, double[] fxs) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(OpticalProperties op, double[] fxs, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(OpticalProperties[] ops, double[] fxs, double[] times) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(OpticalProperties op, double[] fxs, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(OpticalProperties[] ops, double[] fxs, double[] fts) + { + throw new NotImplementedException(); + } + + public double[] ROfRho(OpticalProperties[] ops, double rho) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(OpticalProperties[] ops, double rho, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(OpticalProperties[] ops, double[] rhos, double time) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(OpticalProperties op, double rho, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(OpticalProperties op, double[] rhos, double time) + { + throw new NotImplementedException(); + } + + public double[] ROfRhoAndTime(OpticalProperties[] ops, double rho, double time) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(OpticalProperties[] ops, double rho, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(OpticalProperties[] ops, double[] rhos, double ft) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(OpticalProperties op, double rho, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(OpticalProperties op, double[] rhos, double ft) + { + throw new NotImplementedException(); + } + + public Complex[] ROfRhoAndFt(OpticalProperties[] ops, double rho, double ft) + { + throw new NotImplementedException(); + } + + public double[] ROfFx(OpticalProperties[] ops, double fx) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(OpticalProperties[] ops, double fx, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(OpticalProperties[] ops, double[] fxs, double time) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(OpticalProperties op, double fx, double[] times) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(OpticalProperties op, double[] fxs, double time) + { + throw new NotImplementedException(); + } + + public double[] ROfFxAndTime(OpticalProperties[] ops, double fx, double time) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(OpticalProperties[] ops, double fx, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(OpticalProperties[] ops, double[] fxs, double ft) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(OpticalProperties op, double fx, double[] fts) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(OpticalProperties op, double[] fxs, double ft) + { + throw new NotImplementedException(); + } + + public Complex[] ROfFxAndFt(OpticalProperties[] ops, double fx, double ft) + { + throw new NotImplementedException(); + } + + public double[] FluenceOfRhoAndZ(OpticalProperties[] ops, double[] rhos, double[] zs) + { + throw new NotImplementedException(); + } + + public double[] FluenceOfRhoAndZAndTime(OpticalProperties[] ops, double[] rhos, double[] zs, double[] times) + { + throw new NotImplementedException(); + } + + public Complex[] FluenceOfRhoAndZAndFt(OpticalProperties[] ops, double[] rhos, double[] zs, double[] fts) + { + throw new NotImplementedException(); + } + + public double[] FluenceOfFxAndZ(OpticalProperties[] ops, double[] fxs, double[] zs) + { + throw new NotImplementedException(); + } + + public double[] FluenceOfFxAndZAndTime(OpticalProperties[] ops, double[] fxs, double[] zs, double[] times) + { + throw new NotImplementedException(); + } + + public Complex[] FluenceOfFxAndZAndFt(OpticalProperties[] ops, double[] fx, double[] zs, double[] fts) + { + throw new NotImplementedException(); + } + + + #endregion + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo08UnitTestComparison.cs b/src/Vts.Scripting/MonteCarlo/Demo08UnitTestComparison.cs new file mode 100644 index 000000000..9eb23f83b --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo08UnitTestComparison.cs @@ -0,0 +1,92 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate comparing to well-characterized unit test results +/// +internal class Demo08UnitTestComparison : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 08: run a Monte Carlo simulation and verify results with those in unit tests in Visual Studio + // Spell out all input to ensure same settings as in unit test + + // specify all detector inputs + var rhoRange = new DoubleRange(start: 0, stop: 10, number: 101); + var omegaRange = new DoubleRange(start: 0.05, stop: 1, number: 101); + var timeRange = new DoubleRange(start: 0, stop: 1, number: 101); + var rOfRhoDetectorInput = new ROfRhoDetectorInput { Rho = rhoRange, Name = "ROfRho" }; + var rOfRhoAndOmegaDetectorInput = new ROfRhoAndOmegaDetectorInput { Rho = rhoRange, Omega = omegaRange, Name = "ROfRhoAndOmega" }; + var rOfRhoAndTimeDetectorInput = new ROfRhoAndTimeDetectorInput { Rho = rhoRange, Time = timeRange, Name = "ROfRhoAndTime" }; + + // create a SimulationInput object to define the simulation + var simulationInput = new SimulationInput + { + // specify the number of photons + N = 100, + + OutputName = "results", + + // define a point source beam normally incident at the origin + SourceInput = new DirectionalPointSourceInput + { + SourceType = "DirectionalPoint", + PointLocation = new(x: 0, y: 0, z: 0), + Direction = new(ux: 0, uy: 0, uz: 1), + InitialTissueRegionIndex = 1 + }, + + // define a semi-infinite slab tissue geometry with air-tissue boundary (a bottom air layer is necessary) + TissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 20), // tissue "z" range (slab, 20mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(20, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }, + + // specify all three detectors to run in the single simulation + DetectorInputs = new IDetectorInput[] { rOfRhoDetectorInput, rOfRhoAndOmegaDetectorInput, rOfRhoAndTimeDetectorInput}, + + // specify all simulation options + Options = new SimulationOptions + { + Seed = 0, // -1 will generate a random seed + AbsorptionWeightingType = AbsorptionWeightingType.Analog, + PhaseFunctionType = PhaseFunctionType.HenyeyGreenstein + } + }; + + // create and run the simulation + var simulation = new MonteCarloSimulation(simulationInput); + var simulationOutput = simulation.Run(); + + // plot and compare the results using Plotly.NET + var rOfRhoResults = (ROfRhoDetector)simulationOutput.ResultsDictionary[rOfRhoDetectorInput.Name]; + var rOfRhoAndOmegaResults = (ROfRhoAndOmegaDetector)simulationOutput.ResultsDictionary[rOfRhoAndOmegaDetectorInput.Name]; + var rOfRhoAndTimeResults = (ROfRhoAndTimeDetector)simulationOutput.ResultsDictionary[rOfRhoAndTimeDetectorInput.Name]; + + var unitTestsPass = + // currently 0.95492965855137191 (good) + Math.Abs(rOfRhoResults.Mean[0] - 0.95492965855) < 0.00000000001 && + // currently [0.9548488285246218, -0.008172650857889858] (bad) + Math.Abs((rOfRhoAndOmegaResults.Mean[0, 0] - new Complex(0.95492885, -0.000817329)).Magnitude) < 0.00000001 && + // currently 95.4929658551372 (good) + Math.Abs(rOfRhoAndTimeResults.Mean[0, 0] - 95.492965855) < 0.000000001; + + if (!unitTestsPass) + { + throw new Exception("Unit tests failed"); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo09TransmittanceTallies.cs b/src/Vts.Scripting/MonteCarlo/Demo09TransmittanceTallies.cs new file mode 100644 index 000000000..f30ba9b5b --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo09TransmittanceTallies.cs @@ -0,0 +1,69 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate performing transmittance tallies in Monte Carlo simulations +/// +internal class Demo09TransmittanceTallies : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 09: run a Monte Carlo simulation for transmittance tallies with 1000 photons + + // specify all detector inputs + var rhoRange = new DoubleRange(start: 0, stop: 10, number: 101); + var fxRange = new DoubleRange(start: 0, stop: 0.2, number: 51); + var tOfRhoDetectorInput = new TOfRhoDetectorInput { Rho = rhoRange, Name = "TOfRho" }; + var tOfFxDetectorInput = new TOfFxDetectorInput { Fx = fxRange, Name = "TOfFx" }; + + // create a SimulationInput object to define the simulation + var simulationInput = new SimulationInput + { + // specify the number of photons + N = 1000, + + // define a semi-infinite slab tissue geometry with air-tissue boundary (a bottom air layer is necessary) + TissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 2), // tissue "z" range (slab, 2mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(2, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }, + + // specify all three detectors to run in the single simulation + DetectorInputs = new IDetectorInput[] { tOfRhoDetectorInput, tOfFxDetectorInput} + }; + + // create and run the simulation + var simulation = new MonteCarloSimulation(simulationInput); + var simulationOutput = simulation.Run(); + + // plot the T(rho) results using Plotly.NET + var tOfRhoResults = (TOfRhoDetector)simulationOutput.ResultsDictionary[tOfRhoDetectorInput.Name]; + var tOfRhoLogTransmittance = tOfRhoResults.Mean.Select(r => Math.Log(r)).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (rhoRange.GetMidpoints(), "ρ [mm]", "log(T(ρ)) [mm-2]"); + var chart1 = LineChart(detectorMidpoints, tOfRhoLogTransmittance, xLabel, yLabel, title: "log(T(ρ)) [mm-2]"); + + // plot the T(rho) results using Plotly.NET + var tOfFxResults = (TOfFxDetector)simulationOutput.ResultsDictionary[tOfFxDetectorInput.Name]; + var tOfFxTransmittance = tOfFxResults.Mean.Select(t => t.Magnitude).ToArray(); + var chart2 = LineChart(fxRange.AsEnumerable().ToArray(), tOfFxTransmittance, xLabel, yLabel, title: "T(fx)) [unitless]"); + + if (showPlots) + { + chart1.Show(); + chart2.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/MonteCarlo/Demo10ROfFx.cs b/src/Vts.Scripting/MonteCarlo/Demo10ROfFx.cs new file mode 100644 index 000000000..71bca986e --- /dev/null +++ b/src/Vts.Scripting/MonteCarlo/Demo10ROfFx.cs @@ -0,0 +1,45 @@ +namespace Vts.Scripting.MonteCarlo; + +/// +/// Class using the Vts.dll library to demonstrate performing a Monte Carlo simulation that tallies in the native spatial frequency domain +/// +internal class Demo10ROfFx : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Example 10: run R(fx) detector results + + // create a SimulationInput object to define the simulation + var fxRange = new DoubleRange(start: 0, stop: .5, number: 51); + var detectorInput = new ROfFxDetectorInput { Fx = fxRange, TallySecondMoment = true, Name = "ROfFx" }; + var simulationInput = new SimulationInput + { + // specify the number of photons to run + N = 1000, + + // define a single R(fx) detector at spatial frequencies fx + DetectorInputs = new IDetectorInput[] { detectorInput } + }; + + // create the simulation + var simulation = new MonteCarloSimulation(simulationInput); + + // run the simulation + var simulationOutput = simulation.Run(); + + // plot the results with Plotly.NET + var detectorResults = (ROfFxDetector)simulationOutput.ResultsDictionary[detectorInput.Name]; + var complexReflectance = detectorResults.Mean; + var reflectanceMagnitude = complexReflectance.Select(r => r.Magnitude).ToArray(); + var (detectorMidpoints, xLabel, yLabel) = (fxRange.AsEnumerable().ToArray(), "fx [mm-1]", "R(fx) [unitless]"); + var chart = LineChart(detectorMidpoints, reflectanceMagnitude, xLabel, yLabel, title: "R vs fx [unitless]"); + + if (showPlots) + { + chart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/Program.cs b/src/Vts.Scripting/Program.cs new file mode 100644 index 000000000..7519e8d51 --- /dev/null +++ b/src/Vts.Scripting/Program.cs @@ -0,0 +1,29 @@ +// This is a top-level program that lets you run any of the demos in the Vts.Scripts project. +// It is a good place to start if you want to run a demo, or if you want to write your own code. +using Vts.Scripting; + +// to run a single demo, find the one you want by browsing to them, and run it: +Vts.Scripting.ShortCourse.Demo02AnalogVsContinuousWithReflectance.RunDemo(); +// or: Vts.Scripting.MonteCarlo.Demo01ROfRhoSimple.RunDemo(); +// or: Vts.Scripting.ForwardSolvers.Demo01ROfRhoAndFtSingle.RunDemo(); +// or: Vts.Scripting.ShortCourse.Demo01APhotonCountWithFluence.RunDemo(); + +// or, uncomment one of the following lines to run all demos in a category: +// BatchDemoRunner.RunAllMonteCarloDemos(); +// BatchDemoRunner.RunAllForwardSolverDemos(); +// BatchDemoRunner.RunAllShortCourseDemos(); + +// you can also put any of your own code here, if you want, e.g. commenting out this code will run a Monte Carlo simulation: +// var rhos = new DoubleRange(0.0, 10.0, 101); +// var simulationResults = new MonteCarloSimulation( +// new SimulationInput +// { +// N = 10_000, +// DetectorInputs = new IDetectorInput[] { new ROfRhoDetectorInput { Rho = rhos, Name = "ROfRho" } } +// }).Run(); +// var rhoMidpoints = rhos.GetMidpoints(); +// var detectorResults = (ROfRhoDetector)simulationResults.ResultsDictionary["ROfRho"]; +// LineChart( +// xValues: rhoMidpoints, +// yValues: detectorResults.Mean.Select(r => Math.Log(r)).ToArray(), +// xLabel: "ρ [mm]", yLabel: "log(R(ρ) [mm-2])", title: "log(R(ρ))").Show(); \ No newline at end of file diff --git a/src/Vts.Scripting/ScriptHelper.cs b/src/Vts.Scripting/ScriptHelper.cs new file mode 100644 index 000000000..f088eb2ed --- /dev/null +++ b/src/Vts.Scripting/ScriptHelper.cs @@ -0,0 +1,132 @@ +using Plotly.NET.LayoutObjects; + +namespace Vts.Scripting; + +public static class ScriptHelper +{ + /// + /// Helper extension method that returns an array of midpoints, located halfway between the endpoints of the specified range + /// + /// The range of endpoints + /// The corresponding midpoint outputs + public static double[] GetMidpoints(this DoubleRange endpointRange) + { + var endpoints = endpointRange.AsEnumerable().ToArray(); + if (endpoints.Length < 2) + { + return Array.Empty(); + } + + var midpoints = new double[endpoints.Length - 1]; + for (int i = 0; i < midpoints.Length; i++) + { + midpoints[i] = (endpoints[i + 1] + endpoints[i]) / 2; + } + return midpoints; + } + + /// + /// Helper extension method that returns every nth element of the enumerable, starting at the specified skip index + /// + /// + /// the values being filtered + /// number of values to jump forward at a time + /// number of values to initially skip + /// + public static IEnumerable TakeEveryNth(this IEnumerable values, int n, int skip = 0) => + values.Where((_, i) => (i - skip) % n == 0); + + /// + /// Method to create a standard scatter chart from the specified x and y values using Plotly.NET + /// + /// The values for the x axis. + /// The values for the y axis. + /// The label for the x axis. Optional. + /// The label for the y axis. Optional. + /// The title of the chart. Optional. + /// A `GenericChart` instance representing the scatter chart. + public static Plotly.NET.GenericChart.GenericChart ScatterChart(double[] xValues, double[] yValues, string xLabel = "", string yLabel = "", string title = "") + { + return Chart.Point(xValues, yValues).WithStandardStyling(xLabel, yLabel, title); + } + + /// + /// Method to create a standard scatter chart from the specified x and y values using Plotly.NET + /// + /// The values for the x axis. + /// The values for the y axis. + /// The label for the x axis. Optional. + /// The label for the y axis. Optional. + /// The title of the chart. Optional. + /// A `GenericChart` instance representing the line chart. + public static Plotly.NET.GenericChart.GenericChart LineChart(double[] xValues, double[] yValues, string xLabel = "", string yLabel = "", string title = "") + { + return Chart.Line(xValues, yValues).WithStandardStyling(xLabel, yLabel, title); + } + + /// + /// Fluent helper method to apply standard styling to a chart + /// + /// The `GenericChart` instance to apply styling to. + /// The label for the x-axis. Optional. + /// The label for the y-axis. Optional. + /// The title of the chart. Optional. + /// A `GenericChart` instance with standard styling applied. + private static Plotly.NET.GenericChart.GenericChart WithStandardStyling( + this Plotly.NET.GenericChart.GenericChart chart, string xLabel = "", string yLabel = "", string title = "") + { + // uses Plotly.NET.CSharp.ChartExtensions (adding Plotly.NET to the using statements above will break this) + return chart + .WithTraceInfo(title, ShowLegend: !string.IsNullOrWhiteSpace(title)) + .WithXAxisStyle(Title: Plotly.NET.Title.init(xLabel)) + .WithYAxisStyle(Title: Plotly.NET.Title.init(yLabel)) + .WithLegendStyle(X: 0, Y: 150); + } + + /// + /// Helper method to format a heatmap chart using Plotly.NET. + /// + /// An `IEnumerable` of double arrays specifying the z values of the heatmap. + /// An array of double values specifying the x values of the heatmap. + /// An array of double values specifying the y values of the heatmap. + /// An optional label for the x-axis. Default is an empty string. + /// An optional label for the y-axis. Default is an empty string. + /// An optional title for the chart. Default is an empty string. + /// A `GenericChart` instance representing the heatmap chart. + public static Plotly.NET.GenericChart.GenericChart Heatmap( + IEnumerable values, + double[] x, + double[] y, + string xLabel = "", + string yLabel = "", + string title = "") + { + // attn developers: for reference, the following are the type parameters used in the call to Chart2D.Chart.Heatmap: + // Chart2D.Chart.Heatmap(...) + var chart = Plotly.NET.Chart2D.Chart.Heatmap, double, double, double, string>( + zData: values, + X: x, Y: y, + ReverseScale: false, ReverseYAxis: true, + Transpose: true, + Text: title, + ColorScale: Plotly.NET.StyleParam.Colorscale.Hot + ).WithTraceInfo(title, ShowLegend: !string.IsNullOrWhiteSpace(title)) + .WithLegendStyle(X: 0, Y: 150); + + //var chartLayout = Plotly.NET.GenericChart.getLayout(chart); + //var yAxis = Plotly.NET.Layout.getLinearAxisById(Plotly.NET.StyleParam.SubPlotId.NewYAxis(1)).Invoke(chartLayout); + //yAxis.SetValue("scaleanchor", Plotly.NET.StyleParam.LinearAxisId.NewX(1)); + //chart = Plotly.NET.GenericChartExtensions.WithYAxis(chart, yAxis); + + chart = Plotly.NET.GenericChartExtensions + .WithYAxis(chart, LinearAxis.init( + ScaleAnchor: Plotly.NET.StyleParam.LinearAxisId.NewX(1), AxisType: Plotly.NET.StyleParam.AxisType.Linear)) + .WithXAxisStyle(Title: Plotly.NET.Title.init(xLabel), MinMax: new Tuple(x[0], x[^1])) + .WithYAxisStyle(Title: Plotly.NET.Title.init(yLabel), MinMax: new Tuple(y[0], y[^1])); + + chart = Plotly.NET.GenericChartExtensions + .WithColorbar(chart, title: Plotly.NET.Title.init(title)); + + return chart; + } +} diff --git a/src/Vts.Scripting/ShortCourse/Demo01APhotonCountWithFluence.cs b/src/Vts.Scripting/ShortCourse/Demo01APhotonCountWithFluence.cs new file mode 100644 index 000000000..d40082d5e --- /dev/null +++ b/src/Vts.Scripting/ShortCourse/Demo01APhotonCountWithFluence.cs @@ -0,0 +1,124 @@ +namespace Vts.Scripting.ShortCourse; + +/// +/// Class using the Vts.dll library to demonstrate performing a Monte Carlo simulation +/// to estimate fluence versus radial position and depth for increasing photon counts +/// +internal class Demo01APhotonCountWithFluence : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Short Course Monte Carlo Example 01a: run a Monte Carlo simulation to estimate fluence versus + // radial position and depth for increasing photon counts + + // set up the simulation info that will be constant throughout the set of simulations + + // define a point source beam normally incident at the origin + var sourceInput = new DirectionalPointSourceInput + { + SourceType = "DirectionalPoint", + PointLocation = new(x: 0, y: 0, z: 0), + Direction = new(ux: 0, uy: 0, uz: 1), + InitialTissueRegionIndex = 0 + }; + + // define a semi-infinite slab tissue geometry with air-tissue boundary (a bottom air layer is necessary) + var tissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 100), // tissue "z" range ("semi-infinite" slab, 100mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(100, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }; + + // define a single fluence(rho,z) detector by the endpoints of rho and z bins + var rhoRange = new DoubleRange(start: 0, stop: 10, number: 101); + var zRange = new DoubleRange(start: 0, stop: 10, number: 101); + var detectorInputs = new IDetectorInput[] + { + new FluenceOfRhoAndZDetectorInput { Rho = rhoRange, Z = zRange, TallySecondMoment = true, Name = "FluenceOfRhoAndZ" }, // name can be whatever you want + }; + + // specify simulation options + var options = new SimulationOptions + { + Seed = 0, // -1 will generate a random seed + AbsorptionWeightingType = AbsorptionWeightingType.Analog, + PhaseFunctionType = PhaseFunctionType.HenyeyGreenstein + }; + + // define how many different photon counts to simulate + var numPhotonsArray = new[] { 10, 100, 1000, 10000 }; + + // create an array of simulations, one for each different photon count value + var allSimulations = numPhotonsArray.Select(n => + new MonteCarloSimulation( + new SimulationInput( + numberOfPhotons: n, + outputName: "results", + simulationOptions: options, + sourceInput: sourceInput, + tissueInput: tissueInput, + detectorInputs: detectorInputs) + )).ToArray(); + + // run all the simulations in parallel + var allSimulationOutputs = MonteCarloSimulation.RunAll(allSimulations); + + // collect all the results in helpful arrays + var allFluenceDetectors = allSimulationOutputs.Select(detectorResult => + (FluenceOfRhoAndZDetector)detectorResult.ResultsDictionary["FluenceOfRhoAndZ"]).ToArray(); + var allFluenceMeans = allFluenceDetectors.Select(detector => detector.Mean.ToEnumerable().ToArray()).ToArray(); + var allFluenceSecondMoments = allFluenceDetectors.Select(detector => detector.SecondMoment.ToEnumerable().ToArray()).ToArray(); + + // compute the relative error (standard deviation / mean) for each simulation + var allRelativeErrors = numPhotonsArray.Select((numPhotons, npIdx) => + allFluenceMeans[npIdx].Zip(allFluenceSecondMoments[npIdx], (mean, secondMoment) => + Math.Sqrt((secondMoment - mean * mean) / numPhotons) / mean).ToArray()).ToArray(); + + // plot the results using Plotly.NET + var rhos = rhoRange.GetMidpoints(); + var zs = zRange.GetMidpoints(); + var allRhos = rhos.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + var charts = numPhotonsArray.Select((numPhotons, npIdx) => + { + var fluenceRowsToPlot = allFluenceMeans[npIdx] + .Select(f => Math.Log(f)) // take log for visualization purposes (negative infinity/NaN values won't be rendered ) + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var fluenceDataToPlot = fluenceRowsToPlot.Reverse().Concat(fluenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var fluenceMap = Heatmap(values: fluenceDataToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"log(Φ(ρ, z))"); + + var relativeErrorRowsToPlot = allRelativeErrors[npIdx] + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var relativeErrorDataToPlot = relativeErrorRowsToPlot.Reverse().Concat(relativeErrorRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var relativeErrorMap = Heatmap(values: relativeErrorDataToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"error(ρ, z)"); + + var combined = Chart.Grid(new[]{ fluenceMap, relativeErrorMap }, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + return combined; + }); + + if (showPlots) + { + foreach (var chart in charts) + { + chart.Show(); + } + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ShortCourse/Demo01BAnalogVsDiscreteWithFluence.cs b/src/Vts.Scripting/ShortCourse/Demo01BAnalogVsDiscreteWithFluence.cs new file mode 100644 index 000000000..0d38e5d46 --- /dev/null +++ b/src/Vts.Scripting/ShortCourse/Demo01BAnalogVsDiscreteWithFluence.cs @@ -0,0 +1,150 @@ +namespace Vts.Scripting.ShortCourse; + +/// +/// Class using the Vts.dll library to demonstrate performing a Monte Carlo simulation +/// to estimate fluence versus radial position using different absorption weighting methods +/// +internal class Demo01BAnalogVsDiscreteWithFluence : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Short Course Monte Carlo Example 01b: run a Monte Carlo simulation to estimate fluence versus + // radial position and depth using different absorption weighting methods. This example is similar to Example 01a + // but runs 8 total combos: the same 4 different photon counts, now for each of 2 absorption weighting combos + + // set up the simulation info that will be constant throughout the set of simulations + + // define a point source beam normally incident at the origin + var sourceInput = new DirectionalPointSourceInput + { + SourceType = "DirectionalPoint", + PointLocation = new(x: 0, y: 0, z: 0), + Direction = new(ux: 0, uy: 0, uz: 1), + InitialTissueRegionIndex = 0 + }; + + // define a semi-infinite slab tissue geometry with air-tissue boundary (a bottom air layer is necessary) + var tissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 100), // tissue "z" range ("semi-infinite" slab, 100mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(100, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }; + + // define a single fluence(rho,z) detector by the endpoints of rho and z bins + var rhoRange = new DoubleRange(start: 0, stop: 10, number: 101); + var zRange = new DoubleRange(start: 0, stop: 10, number: 101); + var detectorInputs = new IDetectorInput[] + { + new FluenceOfRhoAndZDetectorInput { Rho = rhoRange, Z = zRange, TallySecondMoment = true, Name = "FluenceOfRhoAndZ" }, // name can be whatever you want + }; + + // define how many different photon count and absorption weighting combos to simulate + var weightingTypeArray = new[] { AbsorptionWeightingType.Analog, AbsorptionWeightingType.Discrete }; + var numPhotonsArray = new[] { 10, 100, 1000, 10000 }; + var simulationTuples = weightingTypeArray.SelectMany(wt => numPhotonsArray.Select(np => (numPhotons: np, weightingType: wt))).ToArray(); + // equivalent to: + //var simulationTuples = new[] + //{ + // (numPhotons: 10, weightingType: AbsorptionWeightingType.Analog), + // (numPhotons: 100, weightingType: AbsorptionWeightingType.Analog), + // (numPhotons: 1000, weightingType: AbsorptionWeightingType.Analog), + // (numPhotons: 10000, weightingType: AbsorptionWeightingType.Analog), + // (numPhotons: 10, weightingType: AbsorptionWeightingType.Discrete), + // (numPhotons: 100, weightingType: AbsorptionWeightingType.Discrete), + // (numPhotons: 1000, weightingType: AbsorptionWeightingType.Discrete), + // (numPhotons: 10000, weightingType: AbsorptionWeightingType.Discrete), + //}; + + // create an array of simulations, one for each different photon count value + var allSimulations = simulationTuples.Select(tuple => + new MonteCarloSimulation( + new SimulationInput( + numberOfPhotons: tuple.numPhotons, + outputName: "results", + simulationOptions: new SimulationOptions + { + Seed = 0, // -1 will generate a random seed + AbsorptionWeightingType = tuple.weightingType, + PhaseFunctionType = PhaseFunctionType.HenyeyGreenstein + }, + sourceInput: sourceInput, + tissueInput: tissueInput, + detectorInputs: detectorInputs) + )).ToArray(); + + // run all the simulations in parallel + var allSimulationOutputs = MonteCarloSimulation.RunAll(allSimulations); + + // collect all the results in helpful arrays + var allFluenceDetectors = allSimulationOutputs.Select(detectorResult => + (FluenceOfRhoAndZDetector)detectorResult.ResultsDictionary["FluenceOfRhoAndZ"]).ToArray(); + var allFluenceMeans = allFluenceDetectors.Select(detector => detector.Mean.ToEnumerable().ToArray()).ToArray(); + var allFluenceSecondMoments = allFluenceDetectors.Select(detector => detector.SecondMoment.ToEnumerable().ToArray()).ToArray(); + + // compute the relative error (standard deviation / mean) for each simulation + var allRelativeErrors = simulationTuples.Select((tuple, tupleIdx) => + allFluenceMeans[tupleIdx].Zip(allFluenceSecondMoments[tupleIdx], (mean, secondMoment) => + Math.Sqrt((secondMoment - mean * mean) / tuple.numPhotons) / mean).ToArray()).ToArray(); + + // compute the relative error difference between the analog and discrete absorption weighting simulations + var analogRelativeErrors = allRelativeErrors.Take(4).ToArray(); + var discreteRelativeErrors = allRelativeErrors.Skip(4).ToArray(); + var relativeErrorDifference = analogRelativeErrors.Zip(discreteRelativeErrors, (analog, discrete) => + analog.Zip(discrete, (a, d) => a - d).ToArray()).ToArray(); + + // plot the results using Plotly.NET + var rhos = rhoRange.GetMidpoints(); + var zs = zRange.GetMidpoints(); + var allRhos = rhos.Select(rho => -rho).Reverse().Concat(rhos).ToArray(); // duplicate for -rho to make symmetric + var fluenceAndRelativeErrorCharts = simulationTuples.Select((tuple, tupleIdx) => + { + var fluenceRowsToPlot = allFluenceMeans[tupleIdx] + .Select(f => Math.Log(f)) // take log for visualization purposes (negative infinity/NaN values won't be rendered ) + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var fluenceDataToPlot = fluenceRowsToPlot.Reverse().Concat(fluenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var fluenceMap = Heatmap(values: fluenceDataToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"log(Φ(ρ, z)) - {tuple.weightingType} N={tuple.numPhotons}"); + + var relativeErrorRowsToPlot = allRelativeErrors[tupleIdx] + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var relativeErrorDataToPlot = relativeErrorRowsToPlot.Reverse().Concat(relativeErrorRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var relativeErrorMap = Heatmap(values: relativeErrorDataToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"error(ρ, z) - {tuple.weightingType} N={tuple.numPhotons}"); + + var combined = Chart.Grid(new[]{ fluenceMap, relativeErrorMap }, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + return combined; + }); + var analogVsDiscreteRelativeErrorCharts = numPhotonsArray.Select((numPhotons, npIdx) => + { + var differenceRowsToPlot = relativeErrorDifference[npIdx] + .Chunk(zs.Length) // break the heatmap into rows (inner dimension is zs) + .ToArray(); + var differenceDataToPlot = differenceRowsToPlot.Reverse().Concat(differenceRowsToPlot).ToArray(); // duplicate for -rho to make symmetric + var differenceMap = Heatmap(values: differenceDataToPlot, x: allRhos, y: zs, + xLabel: "ρ [mm]", yLabel: "z [mm]", title: $"Δ-error(ρ, z)"); + return differenceMap; + }); + + if (showPlots) + { + fluenceAndRelativeErrorCharts.ForEach(chart => chart.Show()); + analogVsDiscreteRelativeErrorCharts.ForEach(chart => chart.Show()); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/ShortCourse/Demo02AnalogVsContinuousWithReflectance.cs b/src/Vts.Scripting/ShortCourse/Demo02AnalogVsContinuousWithReflectance.cs new file mode 100644 index 000000000..c9448ca9c --- /dev/null +++ b/src/Vts.Scripting/ShortCourse/Demo02AnalogVsContinuousWithReflectance.cs @@ -0,0 +1,108 @@ +namespace Vts.Scripting.ShortCourse; + +/// +/// Class using the Vts.dll library to demonstrate performing a Monte Carlo simulation +/// to estimate reflectance versus radial position using different absorption weighting methods +/// +internal class Demo02AnalogVsContinuousWithReflectance : IDemoScript +{ + /// + /// Sample script to demonstrate this class' stated purpose + /// + public static void RunDemo(bool showPlots = true) + { + // Short Course Monte Carlo Example 02: run a Monte Carlo simulation to estimate reflectance versus + // radial position and depth using different absorption weighting methods. This example is similar to Example 01a + // but runs 8 total combos: the same 4 different photon counts, now for each of 2 absorption weighting combos + + // set up the simulation info that will be constant throughout the set of simulations + + // define a point source beam normally incident at the origin + var sourceInput = new DirectionalPointSourceInput + { + SourceType = "DirectionalPoint", + PointLocation = new(x: 0, y: 0, z: 0), + Direction = new(ux: 0, uy: 0, uz: 1), + InitialTissueRegionIndex = 0 + }; + + // define a semi-infinite slab tissue geometry with air-tissue boundary (a bottom air layer is necessary) + var tissueInput = new MultiLayerTissueInput + { + Regions = new ITissueRegion[] + { + new LayerTissueRegion( + zRange: new(double.NegativeInfinity, 0), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)), // air optical properties + new LayerTissueRegion( + zRange: new(0, 100), // tissue "z" range ("semi-infinite" slab, 100mm thick) + op: new(mua: 0.01, musp: 1.0, g: 0.8, n: 1.4)), // tissue optical properties + new LayerTissueRegion( + zRange: new(100, double.PositiveInfinity), // air "z" range + op: new(mua: 0.0, musp: 1E-10, g: 1.0, n: 1.0)) // air optical properties + } + }; + + // define a single fluence(rho,z) detector by the endpoints of rho and z bins + var rhoRange = new DoubleRange(start: 0, stop: 10, number: 101); + var detectorInputs = new IDetectorInput[] + { + new ROfRhoDetectorInput { Rho = rhoRange, TallySecondMoment = true, Name = "ROfRho" }, // name can be whatever you want + }; + + // define how many different absorption weighting types to simulate + var absorptionWeightingTypes = new[] { AbsorptionWeightingType.Analog, AbsorptionWeightingType.Continuous }; + + // create an array of simulations, one for each different photon count value + var allSimulations = absorptionWeightingTypes.Select(weightingType => + new MonteCarloSimulation( + new SimulationInput( + numberOfPhotons: 10000, + outputName: "results", + simulationOptions: new SimulationOptions + { + Seed = 0, // -1 will generate a random seed + AbsorptionWeightingType = weightingType, + PhaseFunctionType = PhaseFunctionType.HenyeyGreenstein + }, + sourceInput: sourceInput, + tissueInput: tissueInput, + detectorInputs: detectorInputs) + )).ToArray(); + + // run all the simulations in parallel + var allSimulationOutputs = MonteCarloSimulation.RunAll(allSimulations); + + // collect all the results in helpful arrays + var allReflectanceDetectors = allSimulationOutputs.Select(detectorResult => + (ROfRhoDetector)detectorResult.ResultsDictionary["ROfRho"]).ToArray(); + var allReflectanceMeans = allReflectanceDetectors.Select(detector => detector.Mean).ToArray(); + var allReflectanceSecondMoments = allReflectanceDetectors.Select(detector => detector.SecondMoment).ToArray(); + + // compute the relative error (standard deviation / mean) for each simulation + var allRelativeErrors = absorptionWeightingTypes.Select((awt, awtIdx) => + allReflectanceMeans[awtIdx].Zip(allReflectanceSecondMoments[awtIdx], (mean, secondMoment) => + Math.Sqrt((secondMoment - mean * mean) / allSimulations[awtIdx].Input.N) / mean).ToArray()).ToArray(); + + // compute the relative error difference between the analog and discrete absorption weighting simulations + var analogRelativeErrors = allRelativeErrors[0]; + var continuousRelativeErrors = allRelativeErrors[1]; + var relativeErrorDifference = analogRelativeErrors.Zip(continuousRelativeErrors, (a, c) => a - c).ToArray(); + + // plot the results using Plotly.NET + var rhos = rhoRange.GetMidpoints(); + var (aawIdx, cawIdx) = (0, 1); // (analog, continuous) absorption weighting indices for convenience + var reflectanceChart = Chart.Combine(new[] + { + LineChart(rhos, allReflectanceMeans[aawIdx].Select(f => Math.Log(f)).ToArray(), xLabel: "ρ [mm]", yLabel: $"log(R(ρ) [mm-2])", title: "Analog"), + LineChart(rhos, allReflectanceMeans[cawIdx].Select(f => Math.Log(f)).ToArray(), xLabel: "ρ [mm]", yLabel: $"log(R(ρ) [mm-2])", title: "CAW") + }); + var errorChart = LineChart(rhos, relativeErrorDifference, xLabel: "ρ [mm]", yLabel: $"Δ-error(ρ) (Analog-CAW)"); + var combinedChart = Chart.Grid(new[]{ reflectanceChart, errorChart }, nRows: 2, nCols: 1, Pattern: Plotly.NET.StyleParam.LayoutGridPattern.Coupled); + + if (showPlots) + { + combinedChart.Show(); + } + } +} \ No newline at end of file diff --git a/src/Vts.Scripting/Vts.Scripting.csproj b/src/Vts.Scripting/Vts.Scripting.csproj new file mode 100644 index 000000000..8d0479775 --- /dev/null +++ b/src/Vts.Scripting/Vts.Scripting.csproj @@ -0,0 +1,20 @@ + + + + Exe + net6.0 + enable + enable + true + preview + + + + + + + + + + + diff --git a/src/Vts.sln b/src/Vts.sln index 6a4b1d868..c1a3609eb 100644 --- a/src/Vts.sln +++ b/src/Vts.sln @@ -27,6 +27,10 @@ Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Mgrte.Console", "Vts.Mg EndProject Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.FemModeling", "Vts.FemModeling\Vts.FemModeling.csproj", "{6BA37669-7BA8-4FF2-911A-0E70EECD450B}" EndProject +Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Scripting", "Vts.Scripting\Vts.Scripting.csproj", "{2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}" +EndProject +Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Scripting.Test", "Vts.Scripting.Test\Vts.Scripting.Test.csproj", "{529AE3EE-28AF-4C0F-8D07-DC89833128DA}" +EndProject Global GlobalSection(SolutionConfigurationPlatforms) = preSolution Benchmark|Any CPU = Benchmark|Any CPU @@ -285,6 +289,70 @@ Global {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|Win32.Build.0 = Release|Any CPU {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|x86.ActiveCfg = Release|Any CPU {6BA37669-7BA8-4FF2-911A-0E70EECD450B}.ReleaseWhiteList|x86.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Any CPU.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Any CPU.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Mixed Platforms.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Mixed Platforms.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Win32.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|Win32.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|x86.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Benchmark|x86.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|Any CPU.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|Win32.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|Win32.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|x86.ActiveCfg = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Debug|x86.Build.0 = Debug|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|Any CPU.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|Any CPU.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|Mixed Platforms.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|Win32.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|Win32.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|x86.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.Release|x86.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|Any CPU.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|Any CPU.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|Mixed Platforms.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|Mixed Platforms.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|Win32.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|Win32.Build.0 = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|x86.ActiveCfg = Release|Any CPU + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}.ReleaseWhiteList|x86.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|Any CPU.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|Any CPU.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|Mixed Platforms.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|Mixed Platforms.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|Win32.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|Win32.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|x86.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Benchmark|x86.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|Any CPU.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|Any CPU.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|Win32.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|Win32.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|x86.ActiveCfg = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Debug|x86.Build.0 = Debug|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|Any CPU.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|Any CPU.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|Mixed Platforms.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|Win32.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|Win32.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|x86.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.Release|x86.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|Any CPU.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|Any CPU.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|Mixed Platforms.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|Mixed Platforms.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|Win32.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|Win32.Build.0 = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|x86.ActiveCfg = Release|Any CPU + {529AE3EE-28AF-4C0F-8D07-DC89833128DA}.ReleaseWhiteList|x86.Build.0 = Release|Any CPU EndGlobalSection GlobalSection(SolutionProperties) = preSolution HideSolutionNode = FALSE @@ -299,6 +367,8 @@ Global {3B0676C4-5CF1-4A0B-97E4-0CDF2B035702} = {D6A3DE7A-7E93-4932-AA8F-C612FD97830E} {197229DB-8CEC-48DD-A953-FBC58D1A3EBC} = {D6A3DE7A-7E93-4932-AA8F-C612FD97830E} {6BA37669-7BA8-4FF2-911A-0E70EECD450B} = {E87433B1-D691-4D36-BB6B-3B766192CC8E} + {2F4FEE64-9215-4B6C-8180-3C9C02BE7C07} = {D6A3DE7A-7E93-4932-AA8F-C612FD97830E} + {529AE3EE-28AF-4C0F-8D07-DC89833128DA} = {920000AB-D955-4399-AD00-78523C740E2A} EndGlobalSection GlobalSection(ExtensibilityGlobals) = postSolution SolutionGuid = {BB1B236A-1BE9-476A-A4B3-8C0F51F1FDA7} diff --git a/src/Vts.sln.DotSettings b/src/Vts.sln.DotSettings new file mode 100644 index 000000000..0527e7de9 --- /dev/null +++ b/src/Vts.sln.DotSettings @@ -0,0 +1,16 @@ + + True + True + True + True + True + True + True + True + True + True + True + True + True + True + True \ No newline at end of file diff --git a/src/Vts/IO/CustomBinaryStreamWriterOfT.cs b/src/Vts/IO/CustomBinaryStreamWriterOfT.cs index bea033c8d..14541ae4c 100644 --- a/src/Vts/IO/CustomBinaryStreamWriterOfT.cs +++ b/src/Vts/IO/CustomBinaryStreamWriterOfT.cs @@ -78,7 +78,7 @@ private void OpenStream() { // guard against directory not existing ahead of time var path = Path.GetDirectoryName(_filename); - if(!string.IsNullOrWhiteSpace(path)) + if (!string.IsNullOrWhiteSpace(path)) { Directory.CreateDirectory(path); } diff --git a/src/Vts/MonteCarlo/Factories/DatabaseWriterFactory.cs b/src/Vts/MonteCarlo/Factories/DatabaseWriterFactory.cs index a30c1a737..184fe6115 100644 --- a/src/Vts/MonteCarlo/Factories/DatabaseWriterFactory.cs +++ b/src/Vts/MonteCarlo/Factories/DatabaseWriterFactory.cs @@ -36,7 +36,7 @@ public static IList GetSurfaceVirtualBoundaryDatabaseWrite { return databaseTypes.Select(v => GetSurfaceVirtualBoundaryDatabaseWriter(v, filePath, outputName)).ToList(); - + } /// /// Static method to instantiate correct PhotonDatabaseWriter given a @@ -104,7 +104,7 @@ public static CollisionInfoDatabaseWriter GetCollisionInfoDatabaseWriter( { case DatabaseType.pMCDiffuseReflectance: return new CollisionInfoDatabaseWriter(VirtualBoundaryType.pMCDiffuseReflectance, - Path.Combine(filePath, outputName, "CollisionInfoDatabase"), + Path.Combine(filePath, outputName, "CollisionInfoDatabase"), tissue.Regions.Count); case DatabaseType.pMCDiffuseTransmittance: return new CollisionInfoDatabaseWriter(VirtualBoundaryType.pMCDiffuseTransmittance, @@ -121,4 +121,4 @@ public static CollisionInfoDatabaseWriter GetCollisionInfoDatabaseWriter( } } -} +}