Skip to content

Commit

Permalink
Add MC examples #1 and #10
Browse files Browse the repository at this point in the history
  • Loading branch information
dcuccia committed Jul 28, 2023
1 parent 2d0e3bc commit 59facb7
Show file tree
Hide file tree
Showing 8 changed files with 229 additions and 65 deletions.
22 changes: 22 additions & 0 deletions src/Vts.Scripting.Utilities/Utilities.cs
Original file line number Diff line number Diff line change
@@ -0,0 +1,22 @@
using Vts.Common;

namespace Vts.Scripting.Utilities;

public static class PlottingExtensions
{
public static double[] GetMidpoints(this DoubleRange endpointRange)
{
var endpoints = endpointRange.AsEnumerable().ToArray();
if (endpoints.Length < 2)
{
throw new ArgumentException("Endpoints must have at least two elements");
}

var midpoints = new double[endpoints.Length - 1];
for (int i = 0; i < midpoints.Length; i++)
{
midpoints[i] = endpoints[i + 1] - endpoints[i];
}
return endpoints;
}
}
13 changes: 13 additions & 0 deletions src/Vts.Scripting.Utilities/Vts.Scripting.Utilities.csproj
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
<Project Sdk="Microsoft.NET.Sdk">

<PropertyGroup>
<TargetFramework>net7.0</TargetFramework>
<ImplicitUsings>enable</ImplicitUsings>
<Nullable>enable</Nullable>
</PropertyGroup>

<ItemGroup>
<ProjectReference Include="..\Vts\Vts.csproj" />
</ItemGroup>

</Project>
74 changes: 74 additions & 0 deletions src/Vts.Scripting/MonteCarlo01ROfRho.cs
Original file line number Diff line number Diff line change
@@ -0,0 +1,74 @@
using Vts.Common;
using Vts.MonteCarlo.Detectors;
using Vts.MonteCarlo.Sources;
using Vts.MonteCarlo.Tissues;
using Vts.MonteCarlo;
using Vts;
using Vts.Scripting.Utilities;
using Plotly.NET.CSharp;

// Example 01: run a simple Monte Carlo simulation with 1000 photons

// set number of photons
var numPhotons = 1000;

// 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.9, 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
var detectorInput = new ROfRhoDetectorInput { Rho = new(0, 40, 201), TallySecondMoment = true };

// create a SimulationInput object to define the simulation
var simulationInput = new SimulationInput
{
N = numPhotons,
SourceInput = sourceInput,
TissueInput = tissueInput,
DetectorInputs = new IDetectorInput[] { detectorInput },
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
var detectorMidpoints = detectorInput.Rho.GetMidpoints();
var logReflectance = simulationOutput.R_r.Select(r => Math.Log(r)).ToArray();
Chart.Point<double, double, string>(
x: detectorMidpoints,
y: logReflectance
).WithTraceInfo("log(R) vs rho [mm-2]", ShowLegend: true)
.WithXAxisStyle<double, double, string>(Title: Plotly.NET.Title.init("rho [mm]"))
.WithYAxisStyle<double, double, string>(Title: Plotly.NET.Title.init("log(R(rho)) [mm-2]"))
.Show();
73 changes: 73 additions & 0 deletions src/Vts.Scripting/MonteCarlo10ROfFx.cs
Original file line number Diff line number Diff line change
@@ -0,0 +1,73 @@
using Vts.MonteCarlo.Detectors;
using Vts.MonteCarlo.Sources;
using Vts.MonteCarlo.Tissues;
using Vts.MonteCarlo;
using Vts;
using Vts.Scripting.Utilities;
using Plotly.NET.CSharp;

// Example 10: run R(fx) detector results

// set number of photons
var numPhotons = 1000;

// 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.9, 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(fx) detector at spatial frequencies fx
var detectorInput = new ROfFxDetectorInput { Fx = new(0, 1, 101), TallySecondMoment = true };

// create a SimulationInput object to define the simulation
var simulationInput = new SimulationInput
{
N = numPhotons,
SourceInput = sourceInput,
TissueInput = tissueInput,
DetectorInputs = new IDetectorInput[] { detectorInput },
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
var detectorMidpoints = detectorInput.Fx.GetMidpoints();
var reflectance = simulationOutput.R_fx.Select(r => r.Magnitude).ToArray();
Chart.Point<double, double, string>(
x: detectorMidpoints,
y: reflectance
).WithTraceInfo("R vs fx [unitless]", ShowLegend: true)
.WithXAxisStyle<double, double, string>(Title: Plotly.NET.Title.init("fx [mm-1]"))
.WithYAxisStyle<double, double, string>(Title: Plotly.NET.Title.init("R(fx) [unitless]"))
.Show();
49 changes: 0 additions & 49 deletions src/Vts.Scripting/MonteCarloDemo.cs

This file was deleted.

5 changes: 4 additions & 1 deletion src/Vts.Scripting/Properties/launchSettings.json
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
{
"profiles": {
"MonteCarloDemo.cs": {
"MonteCarlo01ROfRho.cs": {
"commandName": "Project"
},
"MonteCarlo10ROfFx.cs": {
"commandName": "Project"
}
}
Expand Down
18 changes: 5 additions & 13 deletions src/Vts.Scripting/Vts.Scripting.csproj
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
<Project Sdk="Microsoft.NET.Sdk">
<Project Sdk="Microsoft.NET.Sdk">

<PropertyGroup>
<OutputType>Exe</OutputType>
Expand All @@ -8,22 +8,14 @@
</PropertyGroup>

<ItemGroup>
<PackageReference Include="SmallSharp" Version="1.1.7">
<PrivateAssets>all</PrivateAssets>
<IncludeAssets>runtime; build; native; contentfiles; analyzers; buildtransitive</IncludeAssets>
</PackageReference>
<PackageReference Include="Plotly.NET" Version="4.1.0" />
<PackageReference Include="Plotly.NET.CSharp" Version="0.11.1" />
<PackageReference Include="SmallSharp" Version="1.1.7" />
</ItemGroup>

<ItemGroup>
<ProjectReference Include="..\Vts.Scripting.Utilities\Vts.Scripting.Utilities.csproj" />
<ProjectReference Include="..\Vts\Vts.csproj" />
</ItemGroup>

<ItemGroup>
<None Remove="MonteCarlo.cs" />
</ItemGroup>

<ItemGroup>
<None Remove="MonteCarloDemo.cs" />
</ItemGroup>

</Project>
40 changes: 38 additions & 2 deletions src/Vts.sln
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,9 @@ 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("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Vts.Scripting", "Vts.Scripting\Vts.Scripting.csproj", "{2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}"
Project("{9A19103F-16F7-4668-BE54-9A1E7A4F7556}") = "Vts.Scripting", "Vts.Scripting\Vts.Scripting.csproj", "{2F4FEE64-9215-4B6C-8180-3C9C02BE7C07}"
EndProject
Project("{FAE04EC0-301F-11D3-BF4B-00C04F79EFBC}") = "Vts.Scripting.Utilities", "Vts.Scripting.Utilities\Vts.Scripting.Utilities.csproj", "{84A508E0-099D-4247-90DC-1D821D87EAD2}"
EndProject
Global
GlobalSection(SolutionConfigurationPlatforms) = preSolution
Expand Down Expand Up @@ -319,6 +321,38 @@ Global
{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
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|Any CPU.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|Any CPU.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|Mixed Platforms.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|Mixed Platforms.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|Win32.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|Win32.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|x86.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Benchmark|x86.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|Any CPU.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|Any CPU.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|Mixed Platforms.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|Mixed Platforms.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|Win32.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|Win32.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|x86.ActiveCfg = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Debug|x86.Build.0 = Debug|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|Any CPU.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|Any CPU.Build.0 = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|Mixed Platforms.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|Mixed Platforms.Build.0 = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|Win32.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|Win32.Build.0 = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|x86.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.Release|x86.Build.0 = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|Any CPU.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|Any CPU.Build.0 = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|Mixed Platforms.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|Mixed Platforms.Build.0 = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|Win32.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|Win32.Build.0 = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|x86.ActiveCfg = Release|Any CPU
{84A508E0-099D-4247-90DC-1D821D87EAD2}.ReleaseWhiteList|x86.Build.0 = Release|Any CPU
EndGlobalSection
GlobalSection(SolutionProperties) = preSolution
HideSolutionNode = FALSE
Expand All @@ -333,9 +367,11 @@ 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}
{84A508E0-099D-4247-90DC-1D821D87EAD2} = {D6A3DE7A-7E93-4932-AA8F-C612FD97830E}
EndGlobalSection
GlobalSection(ExtensibilityGlobals) = postSolution
EnterpriseLibraryConfigurationToolBinariesPath = packages\Unity.2.1.505.2\lib\NET35
SolutionGuid = {BB1B236A-1BE9-476A-A4B3-8C0F51F1FDA7}
EnterpriseLibraryConfigurationToolBinariesPath = packages\Unity.2.1.505.2\lib\NET35
EndGlobalSection
EndGlobal

0 comments on commit 59facb7

Please sign in to comment.