A question about the function Region.RayIntersectBoundary(...)

[DWin233]

Take class InfiniteCylinderTissueRegion as an example,

VTS/src/Vts/MonteCarlo/Tissues/InfiniteCylinderTissueRegion.cs

Lines 128 to 136 in cdc7b3b

	var oneIn = this.ContainsPosition(p1);
	var twoIn = this.ContainsPosition(p2);
	// check if ray within cylinder
	if ((oneIn || _onBoundary) && twoIn)
	{
	return false;
	}
	_onBoundary = false; // reset flag

The _onBoundary here is actually obtained when calculating whether p2 is included in the region, so when p2 is on the region boundary, the above branch will inevitably return false.
While p1 is outside the region and p2 is on the boundary (i.e. oneIn = false, _onBoundary= true, twoIn = true), the line between p1 and p2 may have 0 or 1 intersection (excluding p2). Why to ignore the case where there is 1 intersection and to return false?

[hayakawa16]

Thank you for your very good question. First, a possible scenario is that _onBoundary is set by the check of p1 on line 128, and not updated by the check of p2 on line 129. So _onBoundary is not set only by check on p2. Second, in the scenario you describe where p1 is outside, p2 set _onBoundary to be true and also set twoIn to be true, then the ray is going from outside the cylinder to just the surface and so no intersection.

As an aside, the outer Monte Carlo engine is calling this code to determine if a pseudo-collision is to be placed at the intersection of a photon ray and the cylinder. This pseudo-collision is needed to 1) pro-rate the length of the remainder of the track due to optical property change going from outside in or inside out, and 2) change direction if refraction/reflection occurs. Say you have the following scenario, an infinite cylinder with center (0,0,2) and radius=1.99 is specified in the example infile infile_infinite_cylinder_AOfXAndYAndZ.txt. If "Seed"=0, each time you run infile, same simulation is executed. Setting a break point at line 132 you'll see how the outside code plays a part in this logic.

1. a point source launches a photon with Direction (0,0,1) and ray length 0.1197598
2. at break point p1=(0,0,0) (oneIn=false), p2=(0,0,0.1188) (twoIn=true) so intersection and a pseudo-collision is placed at (0,0,0.01), where the ray intersects the cylinder
3. next time at break point p1=(0,0,0.01) pseudo-collision (oneIn=true _onBoundary=true), p2=(0,0,0.1188) (twoIn=true and no change to _onBoundary) and the if at line 132 is true -> no intersection
   Now modify infile and set radius infinite cylinder to be 2-0.1197598=1.8802402, right at cylinder edge, then
4. a point source launches photon with Direction (0,0,1) with ray length 0.1197598, this puts p2 right at cylinder edge
5. at break point p1=(0,0,0) (oneIn=false), p2=(0,0,0.1197598) (twoIn=true _onBoundary=true) and if at line 132 is true -> no intersection
   So the code has to accommodate instances when there is a pseudo-collision at the ray - tissue region intersection and instances when the ray starts or ends right at the surface of a tissue region.

Let me know if this helps or you'd like further clarification. I'm sorry for my lengthy description. For me it is easier to see code in action.

[DWin233]

Thanks for your reply.
What I suspect is that the line from p1 on the outside of the cylinder to p2 on the surface of the cylinder has passed through the inside of the cylinder.
According to the example you gave, the first step of a photon is 0.1197598 by default. Based on infile_infinite_cylinder_AOfXAndYAndZ.txt, set R=0.0397598, Z=0.08 for the infinite cylinder.
There should be an intersection point at p2'=(0,0,0.0402402) on the path from p1 to p2. It will still be true at line 132, i.e. pseudo-collision, but in fact the actual collision should occur at p2'.

[hayakawa16]

I made your edits (R=0.0397598, Z=0.08 for the infinite cylinder) to my infile. There is an intersection point at (0,0,0.0402402). I should have mentioned on my last post that RayIntersectBoundary gets called twice, and the second call adjusts pseudo-collision point sometimes. Try this:

1. at break point p1=(0,0,0) (oneIn=false) and p2=(0,0,0.1197598) (twoIn=true _onBoundary=true) so if at 132 true and exits.
2. step through code at exit to see where the call to this method is coming from. It is in SingleInclusionTissue.GetDistanceToBoundary. Since no intersection this code continues on and sees if a "projected" photon hits the boundary. In this call to RayIntersectBoundary the photon track length S is set to 100mm (this works for most tissue like systems). At break point p1=(0,0,0) (oneIn=false) and p2=(0,0,100) (twoIn=false) so if at 132 is false and the distanceToBoundary is calculated and is 0.0402402. This distance is used in MonteCarloSimulation.MoveToBoundaryCheck tissueDistance variable and photon.Move is called with this distance to move the photon.
3. At this point if you set a break point in Photon.Move line 159 or just step to this point, you'll see that DP.Position.Z is updated to 0.0402402.

One thing I do when debugging is that I set a break point in Photon.Move line 170 and just see how the photon collision points are getting updated. Then periodically I check on the variable History,HistoryData and look at each PhotonDataPoint Position values to make sure the code is putting a pseudo-collision at the tissue region boundary intersections.

I hope this helps. The code is complex in certain circumstances. I feel you are in the trenches with it (which is inspiring!) and so would like to help you understand as best I can.

[DWin233]

I did focus too much on the implementation and return value of the function and neglected the call to the function. Thanks.

[hayakawa16]

You're welcome!