diff --git a/src/Vts.Test/MonteCarlo/Detectors/pMCCAWOneLayerDetectorsTests.cs b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCCAWOneLayerDetectorsTests.cs similarity index 60% rename from src/Vts.Test/MonteCarlo/Detectors/pMCCAWOneLayerDetectorsTests.cs rename to src/Vts.Test/MonteCarlo/Detectors/pMCdMCCAWOneLayerDetectorsTests.cs index 99ae5ce49..f67fe8600 100644 --- a/src/Vts.Test/MonteCarlo/Detectors/pMCCAWOneLayerDetectorsTests.cs +++ b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCCAWOneLayerDetectorsTests.cs @@ -1,409 +1,525 @@ -using System; -using System.Collections.Generic; -using NUnit.Framework; -using Vts.Common; -using Vts.IO; -using Vts.MonteCarlo; -using Vts.MonteCarlo.Detectors; -using Vts.MonteCarlo.Helpers; -using Vts.MonteCarlo.Sources; -using Vts.MonteCarlo.Tissues; -using Vts.MonteCarlo.PostProcessing; -using Vts.MonteCarlo.PhotonData; - -namespace Vts.Test.MonteCarlo.Detectors -{ - /// - /// These tests execute perturbation Monte Carlo (pMC) on a continuous absorption weighting (CAW) - /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same results, and - /// 2) the tally results match the linux results given the same seed - /// mersenne twister STANDARD_TEST. The linux results assumes photon passes - /// through specular and deweights photon by specular. This test starts photon - /// inside tissue and then multiplies result by specular deweighting to match - /// linux results. The linux results are generated using the post-processing code in - /// the g_post subdirectory. - /// - [TestFixture] - public class pMCCAWOneLayerDetectorsTests - { - private SimulationInput _referenceInputOneLayerTissue; - private SimulationOutput _referenceOutputOneLayerTissue; - private double _factor; - private pMCDatabase _databaseOneLayerTissue; - - private readonly List _listOfTestGeneratedFiles = new() - { - "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed - "DiffuseReflectanceDatabase.txt", - "CollisionInfoDatabase", - "CollisionInfoDatabase.txt", - "file.txt", // file that captures screen output of MC simulation - }; - - [OneTimeTearDown] - public void Clear_folders_and_files() - { - // make sure databases generated from previous tests are deleted - foreach (var file in _listOfTestGeneratedFiles) - { - FileIO.FileDelete(file); - } - } - /// - /// Define SimulationInput to describe homogeneous one layer tissue - /// - /// - /// - [OneTimeSetUp] - public void Execute_Monte_Carlo() - { - // delete previously generated files - Clear_folders_and_files(); - - var simulationOptions = new SimulationOptions( - 0, - RandomNumberGeneratorType.MersenneTwister, - AbsorptionWeightingType.Continuous, - PhaseFunctionType.HenyeyGreenstein, - new List() { DatabaseType.pMCDiffuseReflectance }, - false, // track statistics - 0.0, // RR threshold -> 0 = no RR performed - 0); - var sourceInput = new DirectionalPointSourceInput( - new Position(0.0, 0.0, 0.0), - new Direction(0.0, 0.0, 1.0), - 1); - var detectorInputs = new List() - { - new ROfRhoDetectorInput() { Rho=new DoubleRange(0.0, 10.0, 101)}, - new ROfRhoAndTimeDetectorInput() { Rho=new DoubleRange(0.0, 10.0, 101), Time=new DoubleRange(0.0, 1.0, 101)} , - new ROfFxDetectorInput() { Fx=new DoubleRange(0.0, 0.5, 11)}, - new ROfFxAndTimeDetectorInput() { Fx=new DoubleRange(0.0, 0.5, 11), Time=new DoubleRange(0.0, 1.0, 101)} - }; - _referenceInputOneLayerTissue = new SimulationInput( - 100, - "", // can't create folder in isolated storage - simulationOptions, - sourceInput, - new MultiLayerTissueInput( - new ITissueRegion[] - { - new LayerTissueRegion( - new DoubleRange(double.NegativeInfinity, 0.0), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), - new LayerTissueRegion( - new DoubleRange(0.0, 20.0), - new OpticalProperties(0.01, 1.0, 0.8, 1.4)), - new LayerTissueRegion( - new DoubleRange(20.0, double.PositiveInfinity), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) - } - ), - detectorInputs); - _factor = 1.0 - Optics.Specular( - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); - _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); - - _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); - } - - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference - /// - [Test] - public void Validate_pMC_CAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoAndTimeDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - Time=new DoubleRange(0.0, 1.0, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - - _referenceOutputOneLayerTissue.R_rt[0, 0]), 0.00000000001); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] * _factor - 92.2411018), 0.0000001); - } - [Test] - public void Validate_pMC_CAW_ROfRho_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 0.00000000001); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.922411018), 0.000000001); - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(rho) - /// - [Test] - public void Validate_pMC_CAW_ROfRho_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 1.013156), 0.000001); - } - /// - /// Test to validate that calling dMC results in not a NaN - /// - [Test] - public void Validate_dMC_CAW_dROfRhodMua_dROfRhodMus_produces_not_NaN_results() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new dMCdROfRhodMuaDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 11), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - }, - new dMCdROfRhodMusDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 11), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); - } - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference - /// - [Test] - public void Validate_pMC_CAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxAndTimeDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - Time=new DoubleRange(0.0, 1.0, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 0.00000000001); - } - [Test] - public void Validate_pMC_CAW_ROfFx_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - // set perturbed ops to reference ops - PerturbedOps=new List() { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Real - - _referenceOutputOneLayerTissue.R_fx[0].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Imaginary - - _referenceOutputOneLayerTissue.R_fx[0].Imaginary), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.483629), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0 - - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(fx) - /// - [Test] - public void Validate_pMC_CAW_ROfFx_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from prior run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.303789), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.036982), 0.000001); - } - /// - /// Test to validate that calling dMC results in not a NaN - /// - [Test] - public void Validate_dMC_CAW_dROfRhodMua_produces_not_NaN_results() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new dMCdROfRhodMuaDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - }, - new dMCdROfRhodMusDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); - } - } -} - +using System; +using System.Collections.Generic; +using NUnit.Framework; +using Vts.Common; +using Vts.IO; +using Vts.MonteCarlo; +using Vts.MonteCarlo.Detectors; +using Vts.MonteCarlo.Helpers; +using Vts.MonteCarlo.Sources; +using Vts.MonteCarlo.Tissues; +using Vts.MonteCarlo.PostProcessing; +using Vts.MonteCarlo.PhotonData; + +namespace Vts.Test.MonteCarlo.Detectors +{ + /// + /// These tests execute perturbation Monte Carlo (pMC) on a continuous absorption weighting (CAW) + /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same results, and + /// 2) the tally results match the linux results given the same seed + /// mersenne twister STANDARD_TEST. The linux results assumes photon passes + /// through specular and deweights photon by specular. This test starts photon + /// inside tissue and then multiplies result by specular deweighting to match + /// linux results. The linux results are generated using the post-processing code in + /// the g_post subdirectory. + /// + [TestFixture] + public class pMCdMCCAWOneLayerDetectorsTests + { + private SimulationInput _referenceInputOneLayerTissue; + private SimulationOutput _referenceOutputOneLayerTissue; + private double _factor; + private pMCDatabase _databaseOneLayerTissue; + + private readonly List _listOfTestGeneratedFiles = new() + { + "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed + "DiffuseReflectanceDatabase.txt", + "CollisionInfoDatabase", + "CollisionInfoDatabase.txt", + "file.txt", // file that captures screen output of MC simulation + }; + + [OneTimeTearDown] + public void Clear_folders_and_files() + { + // make sure databases generated from previous tests are deleted + foreach (var file in _listOfTestGeneratedFiles) + { + FileIO.FileDelete(file); + } + } + + /// + /// Define SimulationInput to describe homogeneous one layer tissue + /// + /// + /// + [OneTimeSetUp] + public void Execute_Monte_Carlo() + { + // delete previously generated files + Clear_folders_and_files(); + + var simulationOptions = new SimulationOptions( + 0, + RandomNumberGeneratorType.MersenneTwister, + AbsorptionWeightingType.Continuous, + PhaseFunctionType.HenyeyGreenstein, + new List { DatabaseType.pMCDiffuseReflectance }, + false, // track statistics + 0.0, // RR threshold -> 0 = no RR performed + 0); + var sourceInput = new DirectionalPointSourceInput( + new Position(0.0, 0.0, 0.0), + new Direction(0.0, 0.0, 1.0), + 1); + var detectorInputs = new List + { + new ROfRhoDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101)}, + new ROfRhoRecessedDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101),ZPlane=-1.0}, + new ROfRhoAndTimeDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101), Time=new DoubleRange(0.0, 1.0, 101)} , + new ROfFxDetectorInput { Fx=new DoubleRange(0.0, 0.5, 11)}, + new ROfFxAndTimeDetectorInput { Fx=new DoubleRange(0.0, 0.5, 11), Time=new DoubleRange(0.0, 1.0, 101)} + }; + _referenceInputOneLayerTissue = new SimulationInput( + 100, + "", // can't create folder in isolated storage + simulationOptions, + sourceInput, + new MultiLayerTissueInput( + new ITissueRegion[] + { + new LayerTissueRegion( + new DoubleRange(double.NegativeInfinity, 0.0), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), + new LayerTissueRegion( + new DoubleRange(0.0, 20.0), + new OpticalProperties(0.01, 1.0, 0.8, 1.4)), + new LayerTissueRegion( + new DoubleRange(20.0, double.PositiveInfinity), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) + } + ), + detectorInputs); + _factor = 1.0 - Optics.Specular( + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); + _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); + + _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(rho,time) and that + /// R(rho,time) recessed to a height of 0 are equal + /// + [Test] + public void Validate_pMC_CAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoAndTimeDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0.0, 1.0, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 } + }, + new pMCROfRhoAndTimeRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0, 1.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] * _factor - 92.2411018), 0.0000001); + Assert.AreEqual(88, postProcessedOutput.pMC_R_rt_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rtr[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + Assert.AreEqual(88, postProcessedOutput.pMC_R_rtr_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values (0 perturbation) + /// determines results equal to reference for R(rho) and R(rho) recessed + /// when recessed height=0 for pMC and dMC results + /// + [Test] + public void Validate_pMC_dMC_CAW_ROfRho_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new pMCROfRhoRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.922411018), 0.000000001); + // validation value based on previous run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 30.0061), 0.0001); + Assert.AreEqual(88, postProcessedOutput.pMC_R_r_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rr[0] - + _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + Assert.AreEqual(88, postProcessedOutput.pMC_R_rr_TallyCount); + + // validate derivatives with respect to mua and mus with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.612979) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.207432) < 1e-6); + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(rho) + /// + [Test] + public void Validate_pMC_dMC_CAW_ROfRho_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 1.013156), 0.000001); + Assert.IsTrue(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 37.0997) < 1e-4); + Assert.AreEqual(88, postProcessedOutput.pMC_R_r_TallyCount); + // validate derivative values with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.609284) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.192882) < 1e-6); + // and 2nd moments + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r2[0] - 19.5494) < 1e-4); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r2[0] - 5.47322) < 1e-5); + } + + /// + /// Test to validate that calling dMC results in not a NaN + /// + [Test] + public void Validate_dMC_CAW_dROfRhodMua_dROfRhodMus_produces_not_NaN_results() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new dMCdROfRhodMuaDetectorInput + { + Rho = new DoubleRange(0.0, 10, 11), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + }, + new dMCdROfRhodMusDetectorInput + { + Rho = new DoubleRange(0.0, 10, 11), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + // validation value obtained from linux run using above input and seeded the same + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference + /// + [Test] + public void Validate_pMC_CAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxAndTimeDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + Time=new DoubleRange(0.0, 1.0, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, + PerturbedRegionsIndices=new List { 1 } + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 1e-10); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 1e-10); + } + + [Test] + public void Validate_pMC_CAW_ROfFx_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Real - + _referenceOutputOneLayerTissue.R_fx[0].Real), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[0].Imaginary - + _referenceOutputOneLayerTissue.R_fx[0].Imaginary), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.483629), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0= + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(fx) + /// + [Test] + public void Validate_pMC_CAW_ROfFx_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from prior run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.303789), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.036982), 0.000001); + } + + /// + /// Test to validate that calling dMC results in not a NaN + /// + [Test] + public void Validate_dMC_CAW_dROfRhodMua_produces_not_NaN_results() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new dMCdROfRhodMuaDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + }, + new dMCdROfRhodMusDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); + } + } +} + diff --git a/src/Vts.Test/MonteCarlo/Detectors/pMCDAWOneLayerDetectorsTests.cs b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs similarity index 59% rename from src/Vts.Test/MonteCarlo/Detectors/pMCDAWOneLayerDetectorsTests.cs rename to src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs index 1fd46a3b2..5b1ccf058 100644 --- a/src/Vts.Test/MonteCarlo/Detectors/pMCDAWOneLayerDetectorsTests.cs +++ b/src/Vts.Test/MonteCarlo/Detectors/pMCdMCDAWOneLayerDetectorsTests.cs @@ -1,411 +1,549 @@ -using System; -using System.Collections.Generic; -using NUnit.Framework; -using Vts.Common; -using Vts.IO; -using Vts.MonteCarlo; -using Vts.MonteCarlo.Detectors; -using Vts.MonteCarlo.Helpers; -using Vts.MonteCarlo.Sources; -using Vts.MonteCarlo.Tissues; -using Vts.MonteCarlo.PostProcessing; -using Vts.MonteCarlo.PhotonData; - -namespace Vts.Test.MonteCarlo.Detectors -{ - /// - /// These tests execute perturbation Monte Carlo (pMC) on a discrete absorption weighting (DAW) - /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same - /// (no perturbation) results, and 2) the tally results match the linux results given the same seed - /// mersenne twister STANDARD_TEST. The linux results assumes photon passes - /// through specular and deweights photon by specular. This test starts photon - /// inside tissue and then multiplies result by specular deweighting to match - /// linux results. The linux results are generated using the post-processing code in - /// the g_post subdirectory. - /// - [TestFixture] - public class pMCDAWOneLayerDetectorsTests - { - private SimulationInput _referenceInputOneLayerTissue; - private SimulationOutput _referenceOutputOneLayerTissue; - private double _factor; - private pMCDatabase _databaseOneLayerTissue; - - readonly List listOfTestFiles = new List() - { - "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed - "DiffuseReflectanceDatabase.txt", - "CollisionInfoDatabase", - "CollisionInfoDatabase.txt", - "file.txt", // file that captures screen output of MC simulation - }; - - [OneTimeTearDown] - public void Clear_folders_and_files() - { - // make sure databases generated from previous tests are deleted - foreach (var file in listOfTestFiles) - { - FileIO.FileDelete(file); - } - } - - /// - /// Define SimulationInput to describe homogeneous one layer case and generate reference database - /// - /// - [OneTimeSetUp] - public void Execute_Monte_Carlo() - { - var simulationOptions = new SimulationOptions( - 0, - RandomNumberGeneratorType.MersenneTwister, - AbsorptionWeightingType.Discrete, - PhaseFunctionType.HenyeyGreenstein, - new List() {DatabaseType.pMCDiffuseReflectance}, - false, // track statistics - 0.0, // RR threshold -> 0 = no RR performed - 0); - var sourceInput = new DirectionalPointSourceInput( - new Position(0.0, 0.0, 0.0), - new Direction(0.0, 0.0, 1.0), - 1); - var detectorInputs = new List() - { - new ROfRhoDetectorInput() {Rho=new DoubleRange(0.0, 10.0, 101)}, - new ROfRhoRecessedDetectorInput() { Rho=new DoubleRange(0.0, 10.0, 101),ZPlane=-1.0}, - new ROfRhoAndTimeDetectorInput() { Rho = new DoubleRange(0.0, 10.0, 101),Time = new DoubleRange(0.0, 1.0, 101)}, - new ROfFxDetectorInput() {Fx=new DoubleRange(0.0, 0.5, 11)}, - new ROfFxAndTimeDetectorInput() { Fx = new DoubleRange(0.0, 0.5, 11),Time = new DoubleRange(0.0, 1.0, 101)} - - }; - _referenceInputOneLayerTissue = new SimulationInput( - 100, - "", // can't create folder in isolated storage - simulationOptions, - sourceInput, - new MultiLayerTissueInput( - new ITissueRegion[] - { - new LayerTissueRegion( - new DoubleRange(double.NegativeInfinity, 0.0), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), - new LayerTissueRegion( - new DoubleRange(0.0, 20.0), - new OpticalProperties(0.01, 1.0, 0.8, 1.4)), - new LayerTissueRegion( - new DoubleRange(20.0, double.PositiveInfinity), - new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) - } - ), - detectorInputs); - _factor = 1.0 - Optics.Specular( - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); - _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); - - _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); - } - - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(rho,time) and that - /// R(rho,time) recessed to a height of 0 are equal - /// - [Test] - public void Validate_pMC_DAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoAndTimeDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - Time=new DoubleRange(0.0, 1.0, 101), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - }, - new pMCROfRhoAndTimeRecessedDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - Time=new DoubleRange(0, 1.0, 101), - ZPlane=0.0, - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - }, - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - - _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0]*_factor - 61.5238307), 0.0000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_rt_TallyCount); - - // validation value obtained from non-pMC non-recessed run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rtr[0, 0] - - _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); - Assert.AreEqual(89, postProcessedOutput.pMC_R_rtr_TallyCount); - } - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(rho) and R(rho) recessed - /// when Height=0, and R(rho,maxdepth) recessed when Height=0 - /// - [Test] - public void Validate_pMC_DAW_ROfRho_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - }, - new pMCROfRhoRecessedDetectorInput() - { - Rho=new DoubleRange(0.0, 10.0, 101), - ZPlane=0.0, - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - }, - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0]*_factor - 0.615238307), 0.000000001); - // validation value based on previous run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 20.022918), 0.000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); - - // validation value obtained from non-pMC non-recessed run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rr[0] - - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); - Assert.AreEqual(89, postProcessedOutput.pMC_R_rr_TallyCount); - - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(rho) - /// - [Test] - public void Validate_pMC_DAW_ROfRho_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfRhoDetectorInput() - { - Rho=new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.7226588), 0.0000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); - } - - /// - /// Test to validate that calling dMC results in not a NaN - /// - [Test] - public void Validate_dMC_DAW_dROfRhodMua_produces_not_NaN_results() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new dMCdROfRhodMuaDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - }, - new dMCdROfRhodMusDetectorInput() - { - Rho = new DoubleRange(0.0, 10, 101), - // set perturbed ops to reference ops - PerturbedOps = new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices = new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from linux run using above input and seeded the same - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); - Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); - Assert.AreEqual(89, postProcessedOutput.dMCdMua_R_r_TallyCount); - Assert.AreEqual(89, postProcessedOutput.dMCdMus_R_r_TallyCount); - } - - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(fx,time) - /// - [Test] - public void Validate_pMC_DAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxAndTimeDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - Time=new DoubleRange(0.0, 1.0, 101), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 } - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - - _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 0.00000000001); - // validation value obtained from prior run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - 6.858014), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - 0.339772), 0.000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_fxt_TallyCount); - } - /// - /// Test to validate that setting mua and mus to the reference values - /// determines results equal to reference for R(fx) - /// - [Test] - public void Validate_pMC_DAW_ROfFx_zero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - PerturbedOps=new List() { // perturbed ops - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP}, - PerturbedRegionsIndices=new List() { 1 }, - TallySecondMoment = true - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from reference non-pMC run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - _referenceOutputOneLayerTissue.R_fx[1].Real), 0.00000000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - _referenceOutputOneLayerTissue.R_fx[1].Imaginary), 0.00000000001); - // validation value based on previous run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.328865), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.083909), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.467357), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0 - Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); - } - /// - /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) - /// determines results equal to linux results for R(fx) - /// - [Test] - public void Validate_pMC_DAW_ROfFx_nonzero_perturbation_one_layer_tissue() - { - var postProcessor = new PhotonDatabasePostProcessor( - VirtualBoundaryType.pMCDiffuseReflectance, - new List() - { - new pMCROfFxDetectorInput() - { - Fx=new DoubleRange(0.0, 0.5, 11), - // set perturbed ops to reference ops - PerturbedOps=new List() - { - _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, - new OpticalProperties( - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, - _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), - _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP - }, - PerturbedRegionsIndices=new List() {1} - } - }, - _databaseOneLayerTissue, - _referenceInputOneLayerTissue); - var postProcessedOutput = postProcessor.Run(); - // validation value obtained from prior run - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.304018), 0.000001); - Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.029895), 0.000001); - Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); - } - - } -} - +using System; +using System.Collections.Generic; +using NUnit.Framework; +using Vts.Common; +using Vts.IO; +using Vts.MonteCarlo; +using Vts.MonteCarlo.Detectors; +using Vts.MonteCarlo.Helpers; +using Vts.MonteCarlo.Sources; +using Vts.MonteCarlo.Tissues; +using Vts.MonteCarlo.PostProcessing; +using Vts.MonteCarlo.PhotonData; + +namespace Vts.Test.MonteCarlo.Detectors +{ + /// + /// These tests execute perturbation Monte Carlo (pMC) on a discrete absorption weighting (DAW) + /// MC simulation with 100 seeded photons and verify that 1) on-the-fly and pMC produces same + /// (no perturbation) results, and 2) the tally results match the linux results given the same seed + /// mersenne twister STANDARD_TEST. The linux results assumes photon passes + /// through specular and deweights photon by specular. This test starts photon + /// inside tissue and then multiplies result by specular deweighting to match + /// linux results. The linux results are generated using the post-processing code in + /// the g_post subdirectory. + /// + [TestFixture] + public class pMCdMCDAWOneLayerDetectorsTests + { + private SimulationInput _referenceInputOneLayerTissue; + private SimulationOutput _referenceOutputOneLayerTissue; + private double _factor; + private pMCDatabase _databaseOneLayerTissue; + + private readonly List _listOfTestFiles = new() + { + "DiffuseReflectanceDatabase", // name has no "test" prefix, it is generated by the code so name fixed + "DiffuseReflectanceDatabase.txt", + "CollisionInfoDatabase", + "CollisionInfoDatabase.txt", + "file.txt", // file that captures screen output of MC simulation + }; + + [OneTimeTearDown] + public void Clear_folders_and_files() + { + // make sure databases generated from previous tests are deleted + foreach (var file in _listOfTestFiles) + { + FileIO.FileDelete(file); + } + } + + /// + /// Define SimulationInput to describe homogeneous one layer case and generate reference database + /// + /// + [OneTimeSetUp] + public void Execute_Monte_Carlo() + { + // delete previously generated files + Clear_folders_and_files(); + + var simulationOptions = new SimulationOptions( + 0, + RandomNumberGeneratorType.MersenneTwister, + AbsorptionWeightingType.Discrete, + PhaseFunctionType.HenyeyGreenstein, + new List {DatabaseType.pMCDiffuseReflectance}, + false, // track statistics + 0.0, // RR threshold -> 0 = no RR performed + 0); + var sourceInput = new DirectionalPointSourceInput( + new Position(0.0, 0.0, 0.0), + new Direction(0.0, 0.0, 1.0), + 1); + var detectorInputs = new List + { + new ROfRhoDetectorInput {Rho=new DoubleRange(0.0, 10.0, 101), TallySecondMoment = true}, + new ROfRhoRecessedDetectorInput { Rho=new DoubleRange(0.0, 10.0, 101),ZPlane=-1.0}, + new ROfRhoAndTimeDetectorInput { Rho = new DoubleRange(0.0, 10.0, 101),Time = new DoubleRange(0.0, 1.0, 101)}, + new ROfFxDetectorInput {Fx=new DoubleRange(0.0, 0.5, 11)}, + new ROfFxAndTimeDetectorInput { Fx = new DoubleRange(0.0, 0.5, 11),Time = new DoubleRange(0.0, 1.0, 101)} + + }; + _referenceInputOneLayerTissue = new SimulationInput( + 100, + "", // can't create folder in isolated storage + simulationOptions, + sourceInput, + new MultiLayerTissueInput( + new ITissueRegion[] + { + new LayerTissueRegion( + new DoubleRange(double.NegativeInfinity, 0.0), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)), + new LayerTissueRegion( + new DoubleRange(0.0, 20.0), + new OpticalProperties(0.01, 1.0, 0.8, 1.4)), + new LayerTissueRegion( + new DoubleRange(20.0, double.PositiveInfinity), + new OpticalProperties(0.0, 1e-10, 1.0, 1.0)) + } + ), + detectorInputs); + _factor = 1.0 - Optics.Specular( + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP.N, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N); + _referenceOutputOneLayerTissue = new MonteCarloSimulation(_referenceInputOneLayerTissue).Run(); + + _databaseOneLayerTissue = pMCDatabase.FromFile("DiffuseReflectanceDatabase", "CollisionInfoDatabase"); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(rho,time) and that + /// R(rho,time) recessed to a height of 0 are equal + /// + [Test] + public void Validate_pMC_DAW_ROfRhoAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoAndTimeDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0.0, 1.0, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 } + }, + new pMCROfRhoAndTimeRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + Time=new DoubleRange(0, 1.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rt[0, 0]*_factor - 61.5238307), 0.0000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_rt_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rtr[0, 0] - + _referenceOutputOneLayerTissue.R_rt[0, 0]), 1e-10); + Assert.AreEqual(89, postProcessedOutput.pMC_R_rtr_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values (0 perturbation) + /// determines results equal to reference for R(rho) and R(rho) recessed + /// when recessed height=0 for pMC and dMC results + /// + [Test] + public void Validate_pMC_dMC_DAW_ROfRho_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new pMCROfRhoRecessedDetectorInput + { + Rho=new DoubleRange(0.0, 10.0, 101), + ZPlane=0.0, + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + PerturbedOps=new List + { // set perturbed ops to reference ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] - _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0]*_factor - 0.615238307), 0.000000001); + // validation value based on previous run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 20.022918), 0.000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); + + // validation value obtained from non-pMC non-recessed run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_rr[0] - + _referenceOutputOneLayerTissue.R_r[0]), 1e-10); + Assert.AreEqual(89, postProcessedOutput.pMC_R_rr_TallyCount); + + // validate derivatives with respect to mua and mus with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.166157) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.213279) < 1e-6); + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(rho) + /// + [Test] + public void Validate_pMC_dMC_DAW_ROfRho_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfRhoDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + new dMCdROfRhodMuaDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + }, + + new dMCdROfRhodMusDetectorInput + { + Rho=new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new ( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1}, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_r[0] * _factor - 0.7226588), 0.0000001); + Assert.IsTrue(Math.Abs(postProcessedOutput.pMC_R_r2[0] - 28.10877) < 1e-4); + Assert.AreEqual(89, postProcessedOutput.pMC_R_r_TallyCount); + // validate derivative values with prior run + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r[0] + 0.187384) < 1e-6); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r[0] - 0.235936) < 1e-6); + // and 2nd moments + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMua_R_r2[0] - 1.80735) < 1e-5); + Assert.IsTrue(Math.Abs(postProcessedOutput.dMCdMus_R_r2[0] - 5.22420) < 1e-5); + } + + /// + /// Test to validate that calling dMC results in not a NaN + /// + [Test] + public void Validate_dMC_DAW_dROfRhodMua_produces_not_NaN_results() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new dMCdROfRhodMuaDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + }, + new dMCdROfRhodMusDetectorInput + { + Rho = new DoubleRange(0.0, 10, 101), + // set perturbed ops to reference ops + PerturbedOps = new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices = new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from linux run using above input and seeded the same + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMua_R_r[0])); + Assert.AreNotEqual(double.NaN, Math.Abs(postProcessedOutput.dMCdMus_R_r[0])); + Assert.AreEqual(89, postProcessedOutput.dMCdMua_R_r_TallyCount); + Assert.AreEqual(89, postProcessedOutput.dMCdMus_R_r_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(fx,time) + /// + [Test] + public void Validate_pMC_DAW_ROfFxAndTime_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxAndTimeDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + Time=new DoubleRange(0.0, 1.0, 101), + PerturbedOps=new List + { // perturbed ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 } + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Real), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - + _referenceOutputOneLayerTissue.R_fxt[1, 0].Imaginary), 0.00000000001); + // validation value obtained from prior run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Real - 6.858014), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fxt[1, 0].Imaginary - 0.339772), 0.000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_fxt_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the reference values + /// determines results equal to reference for R(fx) + /// + [Test] + public void Validate_pMC_DAW_ROfFx_zero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + PerturbedOps=new List + { // perturbed ops + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP, + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List { 1 }, + TallySecondMoment = true + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + + // validation value obtained from reference non-pMC run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - _referenceOutputOneLayerTissue.R_fx[1].Real), 0.00000000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - _referenceOutputOneLayerTissue.R_fx[1].Imaginary), 0.00000000001); + // validation value based on previous run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.328865), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.083909), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Real - 0.467357), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx2[1].Imaginary - 0.0), 0.000001); // imag of 2nd moment is 0 + Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); + } + + /// + /// Test to validate that setting mua and mus to the perturbed values (mua*2, mus*1.1) + /// determines results equal to linux results for R(fx) + /// + [Test] + public void Validate_pMC_DAW_ROfFx_nonzero_perturbation_one_layer_tissue() + { + var postProcessor = new PhotonDatabasePostProcessor( + VirtualBoundaryType.pMCDiffuseReflectance, + new List + { + new pMCROfFxDetectorInput + { + Fx=new DoubleRange(0.0, 0.5, 11), + // set perturbed ops to reference ops + PerturbedOps=new List + { + _referenceInputOneLayerTissue.TissueInput.Regions[0].RegionOP, + new( + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Mua * 2, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.Musp * 1.1, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.G, + _referenceInputOneLayerTissue.TissueInput.Regions[1].RegionOP.N), + _referenceInputOneLayerTissue.TissueInput.Regions[2].RegionOP + }, + PerturbedRegionsIndices=new List {1} + } + }, + _databaseOneLayerTissue, + _referenceInputOneLayerTissue); + var postProcessedOutput = postProcessor.Run(); + // validation value obtained from prior run + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Real - 0.304018), 0.000001); + Assert.Less(Math.Abs(postProcessedOutput.pMC_R_fx[1].Imaginary - 0.029895), 0.000001); + Assert.AreEqual(89, postProcessedOutput.pMC_R_fx_TallyCount); + } + + [Test] + public void Test_Analog_absorption_weighting_type_throws_argument_exception() + { + var test1 = new DMuaDetectorTest(); + Assert.Throws(() => + { + try + { + test1.SetAbsorbAction(AbsorptionWeightingType.Analog); + } + catch (Exception e) + { + Assert.AreEqual("Analog is not allowed with this detector (Parameter 'awt')", e.Message); + throw; + } + }); + var test2 = new DMusDetectorTest(); + Assert.Throws(() => + { + try + { + test2.SetAbsorbAction(AbsorptionWeightingType.Analog); + } + catch (Exception e) + { + Assert.AreEqual("Analog is not allowed with this detector (Parameter 'awt')", e.Message); + throw; + } + }); + } + + /// + /// Expose protected method in a new class that inherits the class under test + /// + public class DMuaDetectorTest : dMCdROfRhodMuaDetector + { + public new void SetAbsorbAction(AbsorptionWeightingType awt) + { + base.SetAbsorbAction(awt); + } + } + + /// + /// Expose protected method in a new class that inherits the class under test + /// + public class DMusDetectorTest : dMCdROfRhodMusDetector + { + public new void SetAbsorbAction(AbsorptionWeightingType awt) + { + base.SetAbsorbAction(awt); + } + } + + } +} + diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs index 73d04588a..04a99d148 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMuaDetector.cs @@ -1,4 +1,4 @@ -using System; + using System; using System.Collections.Generic; using System.Linq; using System.Runtime.Serialization; @@ -138,8 +138,8 @@ public void Initialize(ITissue tissue, Random rng) TallyCount = 0; // if the data arrays are null, create them (only create second moment if TallySecondMoment is true) - Mean = Mean ?? new double[Rho.Count - 1]; - SecondMoment = SecondMoment ?? (TallySecondMoment ? new double[Rho.Count - 1] : null); + Mean ??= new double[Rho.Count - 1]; + SecondMoment ??= (TallySecondMoment ? new double[Rho.Count - 1] : null); // initialize any other necessary class fields here _perturbedOps = PerturbedOps; @@ -155,19 +155,17 @@ public void Initialize(ITissue tissue, Random rng) /// absorption weighting type protected void SetAbsorbAction(AbsorptionWeightingType awt) { - switch (awt) + // AbsorptionWeightingType.Analog cannot have derivatives so exception + _absorbAction = awt switch { - // note: pMC is not applied to analog processing, - // only DAW and CAW - case AbsorptionWeightingType.Continuous: - _absorbAction = AbsorbContinuous; - break; - case AbsorptionWeightingType.Discrete: - default: - _absorbAction = AbsorbDiscrete; - break; - } + // note: pMC is not applied to analog processing, only DAW and CAW + AbsorptionWeightingType.Continuous => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Discrete => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Analog => throw new ArgumentException(@"Analog is not allowed with this detector", nameof(awt)), + _ => throw new ArgumentOutOfRangeException(typeof(AbsorptionWeightingType).ToString()) + }; } + /// /// method to tally to detector /// @@ -190,58 +188,26 @@ public void Tally(Photon photon) SecondMoment[ir] += photon.DP.Weight * weightFactor * photon.DP.Weight * weightFactor; } - private double AbsorbContinuous(IList numberOfCollisions, IList pathLength, IList perturbedOps) - { - var weightFactor = 1.0; - - // NOTE: following code only works for single perturbed region - foreach (var i in _perturbedRegionsIndices) - { - weightFactor *= - -pathLength[i] * // dMua* factor - (Math.Exp(-perturbedOps[i].Mua * pathLength[i]) / Math.Exp(-_referenceOps[i].Mua * pathLength[i])); // mua pert - if (numberOfCollisions[i] > 0) - { - // the following is more numerically stable - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - } - else - { - weightFactor *= Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * pathLength[i]); - } - } - return weightFactor; - } - - private double AbsorbDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) + /// + /// The following method works for both discrete or continuous absorption weighting + /// + /// photon number of collisions in perturbed region + /// photon path length in perturbed region + /// perturbed optical properties of perturbed region + /// derivative with respect to mus + private double AbsorbContinuousOrDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) { - var weightFactor = 1.0; - - // NOTE: following code only works for single perturbed region - foreach (var i in _perturbedRegionsIndices) - { - if (numberOfCollisions[i] > 0) - { - weightFactor *= - -pathLength[i] * // dMua* factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); - } - else // numberOfCollisions[i] in pert region is 0 - { - weightFactor *= - -pathLength[i] * // dMua* factor - Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * - pathLength[i]); - } - } - return weightFactor; + // NOTE: following code only works for single perturbed region because derivative of + // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that + // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation + var i = _perturbedRegionsIndices[0]; + // rearranged to be more numerically stable + return -pathLength[i] * // dMua* factor + Math.Pow( + _perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - _referenceOps[i].Mus - _referenceOps[i].Mua) * + pathLength[i] / numberOfCollisions[i]), + numberOfCollisions[i]); } /// diff --git a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs index 65db23655..8f14c7cd0 100644 --- a/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs +++ b/src/Vts/MonteCarlo/Detectors/dMCdROfRhodMusDetector.cs @@ -1,4 +1,4 @@ -using System; + using System; using System.Collections.Generic; using System.Linq; using System.Runtime.Serialization; @@ -137,10 +137,10 @@ public void Initialize(ITissue tissue, Random rng) TallyCount = 0; // if the data arrays are null, create them (only create second moment if TallySecondMoment is true) - Mean = Mean ?? new double[Rho.Count - 1]; - SecondMoment = SecondMoment ?? (TallySecondMoment ? new double[Rho.Count - 1] : null); + Mean ??= new double[Rho.Count - 1]; + SecondMoment ??= (TallySecondMoment ? new double[Rho.Count - 1] : null); - // intialize any other necessary class fields here + // initialize any other necessary class fields here _perturbedOps = PerturbedOps; _perturbedRegionsIndices = PerturbedRegionsIndices; _referenceOps = tissue.Regions.Select(r => r.RegionOP).ToList(); @@ -154,19 +154,17 @@ public void Initialize(ITissue tissue, Random rng) /// absorption weighting type protected void SetAbsorbAction(AbsorptionWeightingType awt) { - switch (awt) + // AbsorptionWeightingType.Analog cannot have derivatives so exception + _absorbAction = awt switch { - // note: pMC is not applied to analog processing, - // only DAW and CAW - case AbsorptionWeightingType.Continuous: - _absorbAction = AbsorbContinuous; - break; - case AbsorptionWeightingType.Discrete: - default: - _absorbAction = AbsorbDiscrete; - break; - } + // note: pMC is not applied to analog processing, only DAW and CAW + AbsorptionWeightingType.Continuous => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Discrete => AbsorbContinuousOrDiscrete, + AbsorptionWeightingType.Analog => throw new ArgumentException(@"Analog is not allowed with this detector", nameof(awt)), + _ => throw new ArgumentOutOfRangeException(typeof(AbsorptionWeightingType).ToString()) + }; } + /// /// method to tally to detector /// @@ -189,86 +187,28 @@ public void Tally(Photon photon) SecondMoment[ir] += photon.DP.Weight * weightFactor * photon.DP.Weight * weightFactor; } - private double AbsorbContinuous(IList numberOfCollisions, IList pathLength, IList perturbedOps) + /// + /// The following method works for both discrete or continuous absorption weighting + /// + /// photon number of collisions in perturbed region + /// photon path length in perturbed region + /// perturbed optical properties of perturbed region + /// derivative with respect to mua + private double AbsorbContinuousOrDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) { - var weightFactor = 1.0; - - // NOTE: following code only works for single perturbed region - foreach (var i in _perturbedRegionsIndices) - { - weightFactor *= - Math.Exp(-perturbedOps[i].Mus * pathLength[i]) / Math.Exp(-_referenceOps[i].Mus * pathLength[i]); // Mus pert - if (numberOfCollisions[i] - 1 > 0) - { - weightFactor *= - numberOfCollisions[i] * (1.0 / _referenceOps[i].Mus) * // dMus* 1st factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * - pathLength[i] / (numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1) // minus one here - - pathLength[i] * // dMus* 2nd factor - Math.Pow( - _perturbedOps[i].Mus / _referenceOps[i].Mus * - Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * + // NOTE: following code only works for single perturbed region because derivative of + // Radon-Nikodym product needs d(AB)=dA B + A dB and this does not produce that + // Check for only one perturbedRegionIndices specified by user performed in DataStructuresValidation + var i = _perturbedRegionsIndices[0]; + // rearranged to be more numerically stable + return (numberOfCollisions[i] / _perturbedOps[i].Mus - pathLength[i]) * + Math.Pow(_perturbedOps[i].Mus / _referenceOps[i].Mus * + Math.Exp(-(_perturbedOps[i].Mus + _perturbedOps[i].Mua - + _referenceOps[i].Mus - _referenceOps[i].Mua) * pathLength[i] / numberOfCollisions[i]), - numberOfCollisions[i]); // one here - } - else - { - weightFactor *= Math.Exp(-(_perturbedOps[i].Mus - _referenceOps[i].Mus) * pathLength[i]); - } - } - return weightFactor; - } - - private double AbsorbDiscrete(IList numberOfCollisions, IList pathLength, IList perturbedOps) - { - var weightFactor = 1.0; - - // NOTE: following code only works for single perturbed region - foreach (var i in _perturbedRegionsIndices) - { - if (numberOfCollisions[i] - 1 > 0) - { - weightFactor *= - numberOfCollisions[i]*(1.0/_referenceOps[i].Mus)* // dMus* 1st factor - Math.Pow( - _perturbedOps[i].Mus/_referenceOps[i].Mus* - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - - _referenceOps[i].Mus)* - pathLength[i]/(numberOfCollisions[i] - 1)), - numberOfCollisions[i] - 1) // minus one here - - pathLength[i]* // dMus* 2nd factor - Math.Pow( - _perturbedOps[i].Mus/_referenceOps[i].Mus* - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - - _referenceOps[i].Mus)* - pathLength[i]/numberOfCollisions[i]), - numberOfCollisions[i]); // one here - } - else // number of collisions in pert region is 1 - { - weightFactor *= - numberOfCollisions[i] * (1.0 / _referenceOps[i].Mus) * // dMus* 1st factor - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - - _referenceOps[i].Mus) * - pathLength[i]) - - pathLength[i] * // dMus* 2nd factor - (_perturbedOps[i].Mus / _referenceOps[i].Mus) * - Math.Exp( - -(_perturbedOps[i].Mus + _perturbedOps[i].Mus - _referenceOps[i].Mus - - _referenceOps[i].Mus) * - pathLength[i]); - } - } - return weightFactor; + numberOfCollisions[i]); } - /// /// method to normalize detector results after numPhotons launched ///