Feature/122 add ability to use layered concentric infinite cylinders with refractive index mismatch

matlab/post_processing/loadMCResults.m

@@ -22,6 +22,13 @@
             SurfaceFiber.SecondMoment = SurfaceFiber_txt.SecondMoment;
             SurfaceFiber.Stdev = sqrt((SurfaceFiber.SecondMoment - (SurfaceFiber.Mean .* SurfaceFiber.Mean)) / (json.N)); 
             results{di}.SurfaceFiber = SurfaceFiber;
+        case 'SlantedRecessedFiber'
+            SlantedRecessedFiber.Name = detector.Name;
+            SlantedRecessedFiber_txt = readAndParseJson([datadir slash detector.Name '.txt']);
+            SlantedRecessedFiber.Mean = SlantedRecessedFiber_txt.Mean;              
+            SlantedRecessedFiber.SecondMoment = SlantedRecessedFiber_txt.SecondMoment;
+            SlantedRecessedFiber.Stdev = sqrt((SlantedRecessedFiber.SecondMoment - (SlantedRecessedFiber.Mean .* SlantedRecessedFiber.Mean)) / (json.N)); 
+            results{di}.SlantedRecessedFiber = SlantedRecessedFiber;
         case 'RDiffuse'
             RDiffuse.Name = detector.Name;
             RDiffuse_txt = readAndParseJson([datadir slash detector.Name '.txt']);

matlab/post_processing/load_results_script.m

@@ -15,6 +15,7 @@
 outdir = '.';
 
 show.SurfaceFiber =             1;
+show.SlantedRecessedFiber =     1;
 show.RDiffuse =                 1;
 show.ROfRho =                   1;
 show.ROfRhoRecessed =           1;
@@ -94,7 +95,15 @@
             num2str(results{di}.SurfaceFiber.Mean - 3 * results{di}.SurfaceFiber.Stdev) ' - ' ...
             num2str(results{di}.SurfaceFiber.Mean + 3 * results{di}.SurfaceFiber.Stdev)]);
     end
-
+
+    if isfield(results{di}, 'SlantedRecessedFiber') && show.SlantedRecessedFiber
+        disp(['Total reflectance captured by SlantedRecessedFiber detector: ' num2str(results{di}.SlantedRecessedFiber.Mean)]);
+        disp(['Standard Deviation captured by SlantedRecessedFiber detector: ' num2str(results{di}.SlantedRecessedFiber.Stdev)]);
+        disp(['+/- 3sigma by SlantedRecessedFiber detector: ' ...
+            num2str(results{di}.SlantedRecessedFiber.Mean - 3 * results{di}.SlantedRecessedFiber.Stdev) ' - ' ...
+            num2str(results{di}.SlantedRecessedFiber.Mean + 3 * results{di}.SlantedRecessedFiber.Stdev)]);
+    end
+
     if isfield(results{di}, 'RDiffuse') && show.RDiffuse
         disp(['Total reflectance captured by RDiffuse detector: ' num2str(results{di}.RDiffuse.Mean)]);
     end

src/Vts.Test/MonteCarlo/DataStructuresValidation/TissueInputs/MultiConcentricInfiniteCylinderTissueInputValidationTests.cs

@@ -56,6 +56,7 @@ public void validate_underlying_multilayer_tissue_definition()
             var result = SimulationInputValidation.ValidateInput(input);
             Assert.IsFalse(result.IsValid);
         }
+
         /// <summary>
         /// Test to check that infinite cylinders have non-zero axis definitions.
         /// </summary>
@@ -231,6 +232,7 @@ public void validate_infinite_cylinders_are_within_tissue_layer()
             var result = SimulationInputValidation.ValidateInput(input);
             Assert.IsFalse(result.IsValid);
         }
+
         /// <summary>
         /// Test to check that infinite cylinders refractive index matches refractive index of surrounding layer
         /// </summary>
@@ -276,6 +278,7 @@ public void validate_infinite_cylinders_refractive_index_matches_that_of_surroun
             var result = SimulationInputValidation.ValidateInput(input);
             Assert.IsFalse(result.IsValid);
         }
+
         /// <summary>
         /// Test to check two layer tissue with one layer enclosing concentric infinite cylinders works
         /// </summary>

src/Vts.Test/MonteCarlo/DataStructuresValidation/TissueInputs/SingleInfiniteCylinderTissueInputValidationTests.cs

@@ -0,0 +1,207 @@
+using System.Collections.Generic;
+using NUnit.Framework;
+using Vts.Common;
+using Vts.MonteCarlo;
+using Vts.MonteCarlo.Detectors;
+using Vts.MonteCarlo.Sources;
+using Vts.MonteCarlo.Tissues;
+
+namespace Vts.Test.MonteCarlo.DataStructuresValidation.TissueInputs
+{
+    [TestFixture]
+    public class SingleInfiniteCylinderTissueInputValidationTests
+    {        
+        /// <summary>
+        /// Test to check that underlying MultiLayerTissue is good
+        /// </summary>
+        [Test]
+        public void validate_underlying_multilayer_tissue_definition()
+        {
+            var input = new SimulationInput(
+                10,
+                "",
+                new SimulationOptions(),
+                new DirectionalPointSourceInput(),
+                new SingleInfiniteCylinderTissueInput(
+                    new InfiniteCylinderTissueRegion(
+                        new Position(0, 0, 10),
+                        2.0,
+                        new OpticalProperties()),
+                    // define layer tissues that are incorrect
+                    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(100.0, double.PositiveInfinity),
+                            new OpticalProperties(0.0, 1e-10, 1.0, 1.0))
+                    }
+                ),
+                new List<IDetectorInput>() 
+                {
+                    new FluenceOfXAndYAndZDetectorInput()
+                }
+            );
+            var result = SimulationInputValidation.ValidateInput(input);
+            Assert.IsFalse(result.IsValid);
+        }
+
+        /// <summary>
+        /// Test to check that infinite cylinders have non-zero axis definitions.
+        /// </summary>
+        [Test]
+        public void validate_infinite_cylinder_has_nonzero_radii()
+        {
+            var input = new SimulationInput(
+                10,
+                "",
+                new SimulationOptions(),
+                new DirectionalPointSourceInput(),
+                new SingleInfiniteCylinderTissueInput(
+                    // set one infinite cylinder radius to 0.0
+                    new InfiniteCylinderTissueRegion(
+                        new Position(0, 0, 1), 
+                        0.0, // radius=0
+                        new OpticalProperties()),
+                    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))
+                    }
+                ),
+                new List<IDetectorInput>()
+                {
+                    new FluenceOfXAndYAndZDetectorInput()
+                }
+            );
+            var result = SimulationInputValidation.ValidateInput(input);
+            Assert.IsFalse(result.IsValid);
+        }
+
+        /// <summary>
+        /// Test to check that at least one tissue layer is defined
+        /// </summary>
+        [Test]
+        public void validate_at_least_one_tissue_layer_defined()
+        {
+            var input = new SimulationInput(
+                10,
+                "",
+                new SimulationOptions(),
+                new DirectionalPointSourceInput(),
+                new SingleInfiniteCylinderTissueInput(
+                    new InfiniteCylinderTissueRegion(
+                        new Position(0, 0, 1), 
+                        2.0, 
+                        new OpticalProperties()),
+                    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, double.PositiveInfinity),
+                            new OpticalProperties(0.0, 1e-10, 1.0, 1.0))
+                    }
+                ), 
+                new List<IDetectorInput>()
+                {
+                    new FluenceOfXAndYAndZDetectorInput()
+                }
+            );
+            var result = SimulationInputValidation.ValidateInput(input);
+            Assert.IsFalse(result.IsValid);
+        }
+
+        /// <summary>
+        /// Test to check that infinite cylinder is entirely contained within tissue layer
+        /// </summary>
+        [Test]
+        public void validate_infinite_cylinders_are_within_tissue_layer()
+        {
+            var input = new SimulationInput(
+                10,
+                "",
+                new SimulationOptions(),
+                new DirectionalPointSourceInput(),
+                new SingleInfiniteCylinderTissueInput(
+                    // set infinite cylinder radius to go beyond layer
+                    new InfiniteCylinderTissueRegion(
+                        new Position(0, 0, 1), 
+                        15.0, 
+                        new OpticalProperties()), 
+                    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))
+                    }
+                ), 
+                new List<IDetectorInput>() 
+                {
+                    new FluenceOfXAndYAndZDetectorInput() 
+                }
+            );
+            var result = SimulationInputValidation.ValidateInput(input);
+            Assert.IsFalse(result.IsValid);
+        }
+
+        /// <summary>
+        /// Test to check two layer tissue with one layer enclosing concentric infinite cylinders works
+        /// </summary>
+        [Test]
+        public void validate_infinite_cylinder_within_one_layer_of_multilayer_works()
+        {
+            var input = new SimulationInput(
+                10,
+                "",
+                new SimulationOptions(),
+                new DirectionalPointSourceInput(),
+                new SingleInfiniteCylinderTissueInput(
+                    new InfiniteCylinderTissueRegion(
+                        new Position(0, 0, 5), 
+                        2.0,
+                        new OpticalProperties(1.0, 1.0, 0.8, 1.4)),
+                    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, 10.0),
+                            new OpticalProperties(0.01, 1.0, 0.8, 1.4)),
+                        new LayerTissueRegion(
+                            new DoubleRange(10.0, 20.0),
+                            new OpticalProperties(0.1, 1.5, 0.8, 1.3)),
+                        new LayerTissueRegion(
+                            new DoubleRange(20.0, double.PositiveInfinity),
+                            new OpticalProperties(0.0, 1e-10, 1.0, 1.0))
+                    }
+                ),
+                new List<IDetectorInput>() 
+                { 
+                    new FluenceOfXAndYAndZDetectorInput()
+                }
+            );
+            var result = SimulationInputValidation.ValidateInput(input);
+            Assert.IsTrue(result.IsValid);
+        }
+    }
+}

src/Vts.Test/MonteCarlo/Tissues/InfiniteCylinderTissueRegionTests.cs

@@ -37,32 +37,34 @@ public void Validate_infiniteCylinder_properties()
             Assert.AreEqual(0.8, _infiniteCylinderTissueRegion.RegionOP.G);
             Assert.AreEqual(1.4, _infiniteCylinderTissueRegion.RegionOP.N);
         }
+
         /// <summary>
-        /// Validate method RayIntersectBoundary return correct result
+        /// Validate method ContainsPositions return correct Boolean. ContainsPosition is true if inside
+        /// or *on* boundary.
         /// </summary>
         [Test]
-        public void Verify_RayIntersectBoundary_method_returns_correct_result()
+        public void Verify_ContainsPosition_method_returns_correct_result()
         {
-            var photon = new Photon
-            {
-                DP =
-                {
-                    Position = new Position(0, 0, 3),
-                    Direction = new Direction(1, 0, 0)
-                },
-                S = 2.0 // definitely intersect
-            };
-            var result = _infiniteCylinderTissueRegion.RayIntersectBoundary(photon, out var distanceToBoundary);
-            Assert.AreEqual(true, result);
-            Assert.AreEqual(1.0, distanceToBoundary);
-            photon.S = 0.5; // definitely don't intersect
-            result = _infiniteCylinderTissueRegion.RayIntersectBoundary(photon, out distanceToBoundary);
-            Assert.AreEqual(false, result);
-            Assert.AreEqual(double.PositiveInfinity, distanceToBoundary);
-            photon.S = 1.0; // ends right at boundary => both out and no intersection
-            result = _infiniteCylinderTissueRegion.RayIntersectBoundary(photon, out distanceToBoundary);
-            Assert.AreEqual(false, result);
-            Assert.AreEqual(double.PositiveInfinity, distanceToBoundary);
+            var result = _infiniteCylinderTissueRegion.ContainsPosition(new Position(0, 0, 3.0)); // inside
+            Assert.IsTrue(result);
+            result = _infiniteCylinderTissueRegion.ContainsPosition(new Position(0, 0, 2.0)); // on boundary
+            Assert.IsTrue(result);
+        }
+        /// <summary>
+        /// Validate method OnBoundary return correct Boolean.
+        /// </summary>
+        [Test]
+        public void Verify_OnBoundary_method_returns_correct_result()
+        {
+            // OnBoundary returns false if *exactly* on boundary
+            var result = _infiniteCylinderTissueRegion.OnBoundary(new Position(0, 0, 2.0));
+            Assert.IsFalse(result);
+            // but returns true if outside infinite cylinder which doesn't make sense but it is how code is written
+            // and all unit tests (linux included) are based on this wrong return
+            result = _infiniteCylinderTissueRegion.OnBoundary(new Position(0, 0, 0.5));
+            Assert.IsTrue(result);
+            result = _infiniteCylinderTissueRegion.OnBoundary(new Position(0, 0, 2.0));
+            Assert.IsFalse(result);
         }
 
         /// <summary>
@@ -83,11 +85,39 @@ public void Verify_SurfaceNormal_method_returns_correct_result()
             const double x = 0.3; // pick any x value and determine z on cylinder
             var z = Math.Sqrt(_infiniteCylinderTissueRegion.Radius * _infiniteCylinderTissueRegion.Radius - x * x);
             const double y = 1.11; // pick any y value
-            result = _infiniteCylinderTissueRegion.SurfaceNormal(new Position(x, y, z)); 
+            result = _infiniteCylinderTissueRegion.SurfaceNormal(new Position(x, y, z));
             Assert.IsTrue(Math.Abs(result.Ux - 0.145072) < 1e-6);
             Assert.AreEqual(0, result.Uy); // surface normal on infinite cylinder should *always* have Uz=0
             Assert.IsTrue(Math.Abs(result.Uz + 0.989421) < 1e-6);
         }
 
+        /// <summary>
+        /// Validate method RayIntersectBoundary return correct result
+        /// </summary>
+        [Test]
+        public void Verify_RayIntersectBoundary_method_returns_correct_result()
+        {
+            var photon = new Photon
+            {
+                DP =
+                {
+                    Position = new Position(0, 0, 3),
+                    Direction = new Direction(1, 0, 0)
+                },
+                S = 2.0 // definitely intersect
+            };
+            var result = _infiniteCylinderTissueRegion.RayIntersectBoundary(photon, out var distanceToBoundary);
+            Assert.AreEqual(true, result);
+            Assert.AreEqual(1.0, distanceToBoundary);
+            photon.S = 0.5; // definitely don't intersect
+            result = _infiniteCylinderTissueRegion.RayIntersectBoundary(photon, out distanceToBoundary);
+            Assert.AreEqual(false, result);
+            Assert.AreEqual(double.PositiveInfinity, distanceToBoundary);
+            photon.S = 1.0; // ends right at boundary => both out and no intersection
+            result = _infiniteCylinderTissueRegion.RayIntersectBoundary(photon, out distanceToBoundary);
+            Assert.AreEqual(false, result);
+            Assert.AreEqual(double.PositiveInfinity, distanceToBoundary);
+        }
+
     }
 }