Non-flat layered tissues for MCCL

[janakarana]

Is your feature request related to a problem? Please describe.
For some studies or analyses, certain tissues cannot be assumed as flat layers. For example, skin cannot be assumed to be flat if large source-detector separations are considered. Instead of using a mesh-based Monte Carlo model, I want to see if MCCL can have concave or convex surfaces.

Describe the solution you'd like
I would like to see the possibility of using the existing "MultiConcentricInfiniteCylinder" tissue to get convex and concave tissues (See Figure). To do this, the virtual boundaries should be defined at the inner and outer cylindrical surfaces of a hollow cylinder. Such implementation will allow detectors to capture the photons exiting from the cylinder. We can consider that all detectors are normal to the surface.

Describe alternatives you've considered
Defining all layer boundaries by analytical equations.

[hayakawa16]

I've tried to code this up in fits and starts. At first I tried to write code that has a pseudo-"air" Layer (like air in OPs but very small but non-zero mua) region the surrounds the multi-concentric infinite cylinders. This had problems at the initial interface between air and outer cylinder. Then I tried to create a MultiLayerInfiniteCylinderTissue much like a MultiLayerTissue but with our current InfiniteCylinderTissueRegions. Now I'm realizing that I need to create an associated LayerInfiniteCylinderTissueRegion so that it is a "tube" rather than a solid cylinder to make all of the TissueRegion methods work out like RayIntersectBoundary.

I'm floundering a bit. Any comments are welcome.

[dcuccia]

I'm not sure I appreciate what the main technical problem is, but two morsels come to mind:

1. in the past, when we had "internal" air gaps, we'd used zero absorption and very low scattering with g=1. Is that doable here?

2. not sure why we need to do anything special for the outside. Can we just define a normally-incident cylindrical source that's matched to the tissue?

[hayakawa16]

Hi @dcuccia, thanks for your comments! Per your 1: yes! I agree. We had air gaps before. I was initially sure I could get a layer of "air" with the cylinders inside to work. Since the current fiber reflectance detectors could be only be on the reflectance plane, we put the outermost cylinder close to that plane. @janakarana defined very large cylinders and tried to verify code with a MultiLayerTissue of same geometry. This could not be accomplished. Per your 2: the sources are currently defined so that they can be anywhere. I would need to define a new internal fiber detector that could rest on the surface of the outermost cylinder in the air layer.

@janakarana would like refractive index mismatch air above and below the layered cylinders. He has two problems he is analyzing that has a concave downward reflectance geometry, and another that has a concave upward reflectance geometry with a fiber source and several fiber detectors (https://www.mathsisfun.com/calculus/concave-up-down-convex.html).

Are you possibly available at 1pm tomorrow (@lmalenfant cancelled our meeting because she is out this week but you and I could meet)? If not, do you feel that having an "air" LayerTissueRegion that encompasses a series of InfiniteCylinderRegions (with axis along y-axis), with n>1 and mismatch between cylinders, and innermost cylinder being air should work? This would assume I define a new internal fiber detector.

[hayakawa16]

Now that Issue #136 is closed and methods in MultiConcentricInclusionTissue have been fixed, I can work on this issue now.

When I initially approached this problem, I thought of two ways to solve it:
a) define an "air" LayerRegion (with air layer above and below as usual) and concentric cylinder layers with an innermost cylinder of "air", or
b) define a MultiLayerInfiniteCylinder tissue (much like MultiLayerTissue) with air above and below (inside) the layers of mismatched, concentric infinite cylinders layers.

The ease/complication of the two options arises when thinking about how to define the detectors and the virtual boundaries (VBs) to which they are attached (note we have a one-to-one mapping between VBs and detectors).

In a) the new VB could be based on a tissue region surface. There would need to be a reflectance VB on the surface of the outermost cylinder and a transmittance VB on innermost "air" cylinder, but could be general in case other tissue region surfaces would require VBs in the future. This would work for the figure on the left. For the figure on the right, the reflectance VB would be on the innermost cylinder and transmittance would be on the outermost cylinder. So this might be confusing.
In b) there would two new VBs, DiffuseReflectanceInfiniteCylinder on outermost cylinder and DiffuseTransmittanceInfiniteCylinder. In this case the handing of VBs and detectors would be similar to MultiLayerTissue. However, I would have to modify code in MonteCarloSimulation because it currently assumes for all simulations, there is a reflectance and transmittance plane VB.

The new fiber detector associated with the above VBs:
In a) a "direction" of the detector would need to be defined so that the specification of whether it is capturing cylinder outward flow (positive dot product with surface normal) or inward flow (negative dot product with surface normal).
In b) this is a bit easier because the definition of the VB would handle whether flow is in line with detector or not.

The transport:
In both we would start photons inside the topmost cylinder (outer for left figure, inner for right figure) so no loss due to specular between air layer and cylinder (n=1.4).

Any comments about whether a) or b) would be better solution, or possibly a suggested c)?

[hayakawa16]

@dcuccia, do you have any thoughts on this? I'm thinking option a) would be better direction to go, however I'm not sure how best to define the directional VBs. Possibly based on tissue region boundary is not the best because it is limiting and we might need option to detect away from boundary (like RadianceOfRhoAtZ).

[dcuccia]

I think I'd need to dive in a little deeper here. I'm not sure I appreciate the complexities/differences. Without thinking too hard, it seems to me that both drawing scenarios are represented by the the same tissue definition, but with sources and detectors in different locations. The analogy that comes to mind is that we don't have a TransmissionMultiLayeredTissue and ReflectanceMultiLayeredTissue, we just have a MultiLayeredTissue and different locations for sources and detectors. And if you have a detector, you need a corresponding virtual boundary to detect it (and optionally kill it - if this was a "full tube/lumen", you'd want to keep propagating the photon to the other side of the inner wall). Don't VBs have the concept of directionality built-in? E.g. if we wanted to detect flux (i.e. current) in a specific direction, how do we do that today?

[dcuccia]

If it's a cylindrical VB, the VB could be triggered if the vector is pointing "in" or "out"...

[hayakawa16]

Thanks for your response! True we don't have TransmissionMultiLayeredTissue and ReflectanceMultiLayeredTissue, however we have DiffuseReflectanceVirtualBoundary and DiffuseTransmittanceVirtualBoundary which are defined based on z-axis plane and a direction. So yes, directionality is built into VBs.
For this curved, layered geometry, we have been assuming concentric infinite (along y-axis) cylinders with an air innermost cylinder. These are all within a LayerTissueRegion of air. This means that the full cylinder (tube) is part of the system.
Just had an idea: since the detector in both figures is within the air layer (and over all "tissue" system), do think I could use something like our existing VBType.InternalSurface which instantiates the RadianceVB? This RadianceVB was written solely with dosimetry (i.e. at a z-axis plane and passing it in one direction). Instead, I could have a VBType.InfiniteCylinderSurface with a direction based on if dot product with surface normal (assumed to be outward pointed, away from center of curvature).

[dcuccia]

Thinking back, perhaps we could have defined tissue regions using surfaces, so we could just "re-use" those surfaces for VBs. From that spirit, I like your suggestion of the VBType.InfiniteCylinderSurface. and if you have a convention for "outward" directionality, then if it's a negative number, we know it's pointing inward (would definitely track flux in both directions and accumulate a positive or negative number, for generality).

[hayakawa16]

I guess I should define directionality of VBType in name, e.g. VBType.InfiniteCylinderSurfaceInNormalDirection and InfiniteCylinderSurfaceOutNormalDirection (that's kind of confusing). How about exactly how it would be calculated, InfiniteCylinderSurfaceDotNormalPositive and InfiniteCylinderSurfaceDotNormalNegative?

[hayakawa16]

I'm going to give this a try. I can always refactor names when we review code.
Thank you very much @dcuccia! This discussion help me tremendously!

[dcuccia]

That would be detector dependent, so defer to what Janaka needs, but recommend re-using as much as we can from existing detectors...

[hayakawa16]

@janakarana, I am going to use your 131 branch (the 131a branch has method b) above trials).

[janakarana]

@hayakawa16] Please go ahead and do any changes necessary. Thanks

[hayakawa16]

Sounds good @dcuccia about trying to re-use from existing detectors. I think that is why I didn't like option b) above. When trying to code it up, I was creating a bunch of new classes and modifying the MCSimulation class a lot.

[dcuccia]

No problem. In a perfect world, we don't have to add anything to the core simulator - it should "just" determine if a VB is hit, and then tell the corresponding detectors to do what they want with the info.

[hayakawa16]

I have created two VBs for the two situations above, InfiniteCylinderSurfaceDotNormalPositive and InfiniteCylinderSurfaceDotNormalNegative. Now I'm realizing that our detectors have an implied orientation too, e.g. all ROf... and Reflected... are out top surface and likewise all TOf... and Transmitted... are out bottom surface of MultiLayerTissue.
I think I broke protocol when in the past I added SurfaceFiberDetector. And then subsequently SlantedRecessedFiberDetector was added (and I am about to add InfiniteCylinderSurfaceFiberDetector). These were assumed to be on the top surface of a MultiLayerTissue because IsWithinDetectorAperture uses NegativeZAxis as reference to determine if photon is within NA of these detectors. I think these should have their orientation identified in their name and also in how they get processed.

With the two figures above, "Reflectance" in the detector name implies source and detector on same side of tissue, however I don't think we should use "Reflectance" in the name of these fiber detectors. I could name the new detectors similar to the VBs, e.g. InfiniteCylinderSurfaceDotNormalPositiveFiberDetector and InfiniteCylinderSurfaceDotNormalNegativeFiberDetector or possibly InfiniteCylinderSurfaceFiberDotNormalPositiveDetector (move the "Fiber" closer to "Surface"). I'm going to hold off on correcting the existing SurfaceFiberDetector and SlantedRecessedFiberDetector, but think they should be renamed too to indicate DotNormalPositive orientiation.

Any thoughts?

[dcuccia]

I'd hesitate to have different fiber detectors...why not just let the fiber define which way it is pointing? E.g. if there was a sideways detector, we'd just rotate it, not define a new class.

[hayakawa16]

Good feedback, I'll give that a try. Thanks @dcuccia.

[hayakawa16]

After adding a direction axis for the fiber detector, InDirectionOfFiberAxis (not sure about this name, it is meant to indicate the inflow vector direction) and editing other code, I realized that I could greatly simplify things. I don't need to add InfiniteCylinderSurfaceDotNormalPositiveVirtualBoundary and InfiniteCylinderSurfaceDotNormalNegativeVirtualBoundary, I can use the existing RadianceVirtualBoundary I think. This VB was initially defined to handle the dosimetry detector RadianceOfRhoAtZ. I think I can modify this VB to handle other types of internal radiance detection like InfiniteCylinderSurfaceFiberDetector and for both "in" and "out" using the new property InDirectionOfFiberAxis. I would also have to add another property to InfiniteCylinderSurfaceFiberDetector indicating which tissue region this fiber is attached to. If this is added, then the RadianceVB can determine GetDistanceToVirtualBoundary based on the type of internal radiance detection.

I can see some pros/cons to this idea. RadianceVB is now multi-tasking for different internal detectors, however this is much like DiffuseReflectanceVB which handles all ROf... detectors. The difference is that the surface that the detectors "attach" to will be different for different internal radiance detectors specified.

A benefit is that if I can get this to work, further simplification can be done in which we can a single (internal)SurfaceFiberDetector and it can attach to any tissue region surface in either direction.

I'm going to try to get this to work. My prior work with the new VBs was pushed on this branch so I can toss this if it doesn't work.

[hayakawa16]

My idea above (posted 4 days ago) that was inspired by the suggestion to add a Direction of the fiber, seems to be working. I have added unit tests that qualitatively appear to be making sense. Because the RadianceVB can have tallies in either direction of photons at the surface, it is up to the detector to determine, based in the Direction of the fiber, whether to tally or not. Not sure this is the best design, however it means not having duplicate VBs for directional detectors off the same surface.

Not ready for a PR yet, I would like to do some more testing.

[dcuccia]

Sounds like a great design - detectors should always decide if the photon trajectory is within their "aperture".

[hayakawa16]

Thanks for your positive feedback. Makes me feel more confident moving forward.
My only hesitation was that DiffuseReflectanceVB and SpecularReflectionVB share same surface however have different VBs.

[hayakawa16]

@janakarana, I think I'm ready for you to try the 131... branch (not the 131a... branch). In this branch I realized that the unit tests should have the class name that is being tested in the unit test with "Tests" appended. So now there are SurfaceFiberDetectorTests, SlantedRecessedFiberDetectorTests and the newly added InternalFiberDetectorTests that test the new InternalFiberDetector. Please check out the InternalFiberDetectorTests for how I specified a single "tissue" layer with two infinite cylinders and fiber detectors on the top of the outer cylinder and inside the inner cylinder (corresponding to your left and right figures above, respectively).

[hayakawa16]

@janakarana tried putting an internal fiber detector on a layered region. This was not handled correctly by the code so I am in process of trying to correct code. The problem that I'm having is that for infinite cylinders the fibers could be on the outside (dot product with normal >0) or inside (dot product with normal < 0), i.e. a single continuous surface defining region. This isn't the problem, code handled this. However with LayerTissueRegion, there are 2 possibly surfaces that the fiber could be on the inside or outside, the Start layer and the Stop layer. So I'm trying to determine a way to code RadianceVirtualBoundary to handle if the fiber is adjacent to a LayerTissueRegion and then which plane. The detector aperture should handle whether photon is detected (i.e. whether fiber is on inside or outside).

[hayakawa16]

This is an aside question but one that arises when writing unit tests for fiber detectors. SimulationOutput was created just for unit testing, it is not used in the running of MCCL. For unit testing if I would like to test multiple locations of fiber detectors, I have to run a separate simulation for each. This is because to gain access to the detector within the ResultsDictionary we use the linq statement:
public double SurFib { get { return (double)((dynamic)ResultsDictionary[_detectorResults.First(d =>d.TallyType == "SurfaceFiber").Name]).Mean; } }
I would like to specify several fiber detectors within a single simulation, specify each with a unique "Name" and access each result by the "Name". Is that possible? It would speed up unit test time too.
I've tried modifying "double" to "double[]" and variations of linq statements with no luck.

[dcuccia]

Man, cringing at my own decades-old dynamic code (sorry!). Maybe something like this?

public double[] AllSurfaceFiberDetectorMeans
{ 
    get 
    { 
        string[] surfaceFiberDetectorNames = ["UniqueName1", "UniqueName2", "UniqueName3"]; // or can be defined elsewhere
        var surfaceFiberDetectors = surfaceFiberDetectorNames
            .Select(name => { ResultsDictionary.TryGetValue(name, out var sfd); return sfd as SurfaceFiberDetector; })
            .Where(sfd => sfd is not null);
        return surfaceFiberDetectors.Select(sfd => sfd.Mean).ToArray();
    }
}

Or honestly, not sure you need ResultsDictionary for this. The following simpler code would work, too:

public double[] AllSurfaceFiberDetectorMeans2
{ 
    get 
    { 
        string[] surfaceFiberDetectorNames = ["UniqueName1", "UniqueName2", "UniqueName3"]; // or can be defined elsewhere
        IDetector[] surfaceFiberDetectors = _detectorResults
            .Where(detectorResult => surfaceFiberDetectorNames.Contains(detectorResult.Name) && detectorResult is SurfaceFiberDetector)
            .ToArray();
        return surfaceFiberDetectors.Select(sfd => ((SurfaceFiberDetector)sfd).Mean).ToArray();
    }
}

(If we can't support C#12 collection initializers yet, then the string array brackets on the rhs should be replaced by braces.)

[hayakawa16]

Thank you @dcuccia! I will give those a try. You are the king of linq!

[dcuccia]

Haha, thanks

[hayakawa16]

I used the 2nd method described above and it worked quite nicely. Many thanks to the king!

[dcuccia]

Glad it helped :)

[hayakawa16]

I think I'm ready to have this branch reviewed. @janakarana, do you have an infile I could try out before I ask for the review?

[janakarana]

Thanks for doing this. Here are two infiles. One is for flat layers and the other is for concentric cylindrical layers with very large radius.
Expectations:

1. The two results are not significantly different
2. The simulation time of the cylindrical layered case is not significantly longer than that of the flat layered case.

infile_1310_Flat_7Layer_2mmGap_SurfaceDetector.txt
infile_1310_Flat_7Layer_2mmGap_InternalSurfaceDetector.txt

[hayakawa16]

Hi @janakarana, these infiles both specify MultiLayerTissue. Did you attach the correct ones?

Also, one thing I should point out, the dimensions assumed in some of the code takes tissue-like dimensions into account, e.g. slab thicknesses~=100mm. So very large cylinders (e.g. radius~=45,000mm) cannot be handled unless I modify a line of code. I can do this though in a local version of the code and try your infile.

[janakarana]

Hi @hayakawa16
Sorry, I forgot to include this infile.
infile_1310_Cyl_R1_7Layer_2mmGap.txt

I agree that the tissue thickness is unrealistic. We want to examine the effect of curvature on reflection. A very large radius was chosen to mimic the flat layer case. Then, the radius gradually decreases to increase the curvature.

[hayakawa16]

The newly added unit tests that use the nice SimulationOutput additions using the AllDetector specifications run fine on Windows, but fail on Linux. I slightly modified the code above and pulled out the string specification so that it could be used by Means, SecondMoments and TallyCounts. In SimulationOutput like this:

    private readonly string[] _internalSurfaceFiberDetectorNames = { "InternalSurfaceFiber1", "InternalSurfaceFiber2", "InternalSurfaceFiber3", "InternalSurfaceFiber4" };
    /// <summary>
    /// Created array to allow multiple fiber detectors to be specified in single simulation unit test
    /// </summary>
    public double[] AllInternalSurfaceFiberDetectorMeans
    {
        get
        {
           var internalSurfaceFiberDetectors = _detectorResults
                .Where(detectorResult => _internalSurfaceFiberDetectorNames.Contains(detectorResult.Name) &&
                                         detectorResult is InternalSurfaceFiberDetector).ToArray();
            return internalSurfaceFiberDetectors.Select(isfd => ((InternalSurfaceFiberDetector)isfd).Mean).ToArray();
        }
    }
    /// <summary>
    /// Created SecondMoment array to allow multiple fiber detectors to be specified in single simulation unit test
    /// </summary>
    public double[] AllInternalSurfaceFiberDetectorSecondMoments
    {
        get
        {
            var internalSurfaceFiberDetectors = _detectorResults
                .Where(detectorResult => _internalSurfaceFiberDetectorNames.Contains(detectorResult.Name) &&
                                         detectorResult is InternalSurfaceFiberDetector).ToArray();
            return internalSurfaceFiberDetectors.Select(isfd => ((InternalSurfaceFiberDetector)isfd).SecondMoment).ToArray();
        }    /// <summary>
    /// Created TallyCount array to allow multiple fiber detectors to be specified in single simulation unit test
    /// </summary>
    public long[] AllInternalSurfaceFiberDetectorTallyCounts
    {
        get
        {
             var internalSurfaceFiberDetectors = _detectorResults
                .Where(detectorResult => _internalSurfaceFiberDetectorNames.Contains(detectorResult.Name) &&
                                         detectorResult is InternalSurfaceFiberDetector).ToArray();
            return internalSurfaceFiberDetectors.Select(isfd => ((InternalSurfaceFiberDetector)isfd).TallyCount).ToArray();
        }
    }

}
Then on the unit test side, I take note of which fiber is in array location 0,1,2,3 and after I run the simulation, I have code like:
Assert.Less(Math.Abs(_outputCylinderSystem.AllInternalSurfaceFiberDetectorMeans[0] - 0.090464), 0.000001);

Anything pop out that what I'm doing that is wrong?

[hayakawa16]

I did some debugging of this. Turns out that for this unit test, photon #55, track #637, in RayIntersectBoundary, on Windows root1=-0.1666 (Linux agrees), root2=1.13746e-17 (Linux has negative of this number), which means c=-2.84217e-14 (Linux has negative of this number). These inconsistencies means that the photon trajectory on Linux starts to diverge at this track and so the unit test results are different on Linux and do not pass the unit tests (4 of them).

I will try to continue to debug this branch code, however until this is resolved we cannot merge branch into master.

[hayakawa16]

I have returned to trying to figure this out. I found this:
https://learn.microsoft.com/en-us/answers/questions/1097724/net6-c-linux-vs-windows-different-floating-point-c
which provides the answer to the question of Net6 C#, Linux vs Windows, different floating point calculations?

the roslyn compiler (and thus C#) does not implement strict IEEE floating point precision, but rather optimizes for performance. It its not surprising that they don't match when log operations are done. the machines may need to be the
exact same hardware and matching versions of Roslyn (as code optimization can change the precision) to get a match. But
even this may not work.

even though log operations are not involved in the code in question, it makes me feel like I'm never going to get Windows and Linux to agree on this unit test.
Other optiions I have tried or considered:

1. forcing C# to implement strict IEEE floating point rules -> can't find a way to do this yet
2. these small values come about because the quadratic equation is being solved when a ray intersects an infinite cylinder. Rather than the usual solution formula, I found there is a Citardauq (quadratic spelled backwards) Formula (https://en.wikipedia.org/wiki/Quadratic_formula):
   x = (2c)/[-b -/+ sqrt(b*b-4ac)]
   which is an alternative solution formula for quadratic equation. I tried this and unit tests still don't pass, and many others additionally don't pass.
3. The small root2 floating point number, 1.13746e-17, comes about because |4ac| is small so root2 of regular quadratic is b-sqrt(b*b). So on Windows this value is positive and on Linux this value is negative but because in the next set of code I check if (0<root1<1) or (0<root2<1) [which checks if there is an intersection of the ray with the cylinder], then the code takes different paths on each machine. I tried modifying code in MonteCarlo/Helpers/CylinderTissueRegionToolbox/RayIntersectInfiniteCylinder to test if roots are small and negative and if so change their values, but can't get this to work yet.

Any ideas?

[dcuccia]

I was surprised by this. Here's what my chatbot suggested as possibilities for 1., but nothing easy:

https://chatgpt.com/share/672e8db5-677c-8008-8ba6-ea71fadf7e0c

[hayakawa16]

Thank you @dcuccia! It might not be for permanent implementation, but I'll try some of the suggestions. If they work, that would be interesting to know.