The incident angle of remitted photons at tissue surface__NA__ROfRho_MCCL

[future...@gmail.com]

Hello,

I am considering the numerical aperture(NA) of a realistic fiber. When the incident angle of photon that is remitted back to tissue surface is less than collection angle of an optical fiber, it can be collected, or else it is not. 

Then I would like know whether ROfRho of MCCL captures the incident angle of remitted photons at tissue surface. I may need add NA effect by modifying current source code of MCCL.

Thanks.

Owen.

[Carole Hayakawa]

Hi Owen,

Thanks for your question.  We actually have this coded up on a branch of our repository.  The branch is named ch-141202-detector-NA-changes. I am going to consult with our GitHub expert and  either post directions on how to obtain a copy of this branch source code, or a zip of the branch MCCL.

I post again soon.

Best,

Carole

[Carole Hayakawa]

Hi Owen,

I created a Monte Carlo CommandLine (MCCL) zip file using the code from the branch and is too large to attach here so I put it in my google drive:

https://drive.google.com/file/d/1dv92keh04jmAJI8iRs2stwjj7Ouppak8/view?usp=sharing

The version is 0.0.0 because it is an unofficial release of the MCCL.  You'll notice that in the sample infiles, that for the reflectance and transmittance detectors, there is now the ability to specify NA.

Please give it a try and please let me know how it goes for you.

Best,

Carole

[future...@gmail.com]

Hi, Carole,

I did a test run this afternoon. Please see attached .rar file.

The problem is:

the returned   results.ROfRho.Mean is a vector of zeors.

I tried two ways of specifying NA as 0.22, as shown below. But both show the same results, all zeros. 

---------------------------

      "NA": "0.22",

---------------------------

or 

---------------------------

      "NA": 0.22,

---------------------------

Please advise the correct way of specifying NA in infile.

Thanks.

[future...@gmail.com]

in original sample infiles,         "NA": "Infinity",

The definition of NA in source code may be different from the one we commonly use.

On Wednesday, March 6, 2019 at 6:47:50 PM UTC-6, Carole Hayakawa wrote:

[Carole Hayakawa]

Hi Owen,

Using the value without quotes, "NA": 0.22, should work.  I tried your test infile and get the same 0 results.  The unit tests for the NA code all run successfully, so it has to do with running outside of Visual Studio.  I will debug and get back to you.

Thanks for your feedback!

Carole

[future...@gmail.com]

Hi, Carole,

Just share another test finding with you, and hope it is helpful to your debugging.

If NA is set as a value larger than 1.1, the returned results.ROfRho.Mean is same with the results from the case "NA": "Infinity",.

If NA is set as a value less than and equal to 1.0, the returned results.ROfRho.Mean are all Zeros.

I set NA without quotes.

---------------------------

      "NA": 0.22,

---------------------------

[Carole Hayakawa]

Dear Owen,

I found a bug in the ROfRhoDetector code that checks whether the photon is within in NA when the detector FinalTissueRegionIndex is set to 0 (air above the tissue). Only this detector has the bug.  The get-around is to set the FinalTissueRegionIndex to 1 and you should get non-zero reflectance results.  Let me know if this works for you.

Since you found this bug, I'd like to do some more testing before I create a new zip file for you.

Thanks for helping me improve our code!  If you are able to exercise this code and help us flush out any problems/concerns, we plan to merge this branch into master and make it part of our official releases.  So please don't hesitate to give us your comments.

Best,

Carole

[future...@gmail.com]

Carole,

I am glad to test this code for you!

I did a test by setting FinalTissueRegionIndex to 1 in above high-scattering sample (musp = 1000 mm^-1),

And the results show that:

When NA = 0.22,     the total diffuse reflectance collected from a circular area with diameter of 1mm, is 0.0996.

 When NA= Infinite,  the total diffuse reflectance collected from a circular area with diameter of 1mm, is 0.9482

The collection efficiency at this high scattering case =   0.0996/0.9482 = 10%

However, 2.47% is expected for high scattering case,  reported by literature[1].   

2.47% = (NA/n_tiss)^2 = (0.22/1.40)^2,    n_tiss  is refraction index of tissue.

[1] Bargo, P. R., et al. (2003). "Collection efficiency of a single optical fiber in turbid media." Appl Opt 42(16): 3187-3197..

---------------

I replace the directianal point source with a circular Flat source of radius = 0.5mm, the total diffuse reflectance = 0.1068

which results a even higher collection efficiency = 11.26%.

Owen

[Carole Hayakawa]

Dear Owen,

Thank you for the very helpful reference.  After reading it, I ran a test case and got results that don't match those in the paper.  So I returned to the code and added in more unit tests to test the details of the method IsWithinNA.  I found a bug!  I wasn't squaring cos(theta) to obtain sin(theta).  So I fixed that.

I ran a test case from the paper with 1e5 photons, tissue with OPs mua=0.01/mm, musp=1.0/mm, g=0.8, n=1.35, a 0.6 mm detector ring, and two ROfRho detectors, one with NA="Infinity" and one with NA=0.22.  After processing the results, I get an nc=5%.  Now looking at the paper Figure 8(c) [their units are in cm], it looks like they determine nc to be about 3.5%.  I'm not sure how to rectify that difference.  I have found in the past that sometimes it is difficult to replicate results from a paper, either the paper setup is slightly different or some details are missing.  Since I found the two bugs since I create the initial zip, I have generated a new zip file and uploaded it to my google docs.

A couple of comments:

1) In the infile, if you set TallySecondMoment to true, you can bring up the Stdev of the tallies and make sure that Stdev/Mean <= about 5%.  This will ensure you have run enough photons and the statistical errors in your results are not overshadowing the mean results. Running 1e5 in the test case above my Stdev/Mean was 5%.

2) I notice that you are running a very high scattering case (mus'=1000/mm=10000/cm) which is much larger the musp that were executed in this paper (I think max is 40/cm). I wonder if that might cause your results to be different than those in the paper.

Thank you very much for helping to debug this code.  Please let me know if you think the code still isn't working correctly.

Best,

Carole

[future...@gmail.com]

Hi, Carole,

Happy to hear it!  The number returned by your test makes sense to me. I will bring more feedback here after my test on your new code.

Could you send me a link to download the new Zip file you created?  

I downloaded a Zip file from the google drive link you shared in previous post, but after I unzipped it, it looks the  same code as before.  The "Date modified" of  the mc.exe is still 3/6/2019 6:21pm.  All other files also have the same "Date modified".

Cheers!

[Carole Hayakawa]

Hi Owen,

Sorry about that!  Here is the link:

https://drive.google.com/file/d/1U-M47c3LhW-PHqJ3jFeGLKBiPGEIbTvG/view?usp=sharing

Thanks for your help!

Carole

[future...@gmail.com]

Hi, Carole,

 

the attached figure is my calculated light collection efficiency due to limited numerical aperture.

The attached musp_p1_NA.rar is the spatial result when musp = 0.1 mm^-1; 

Compared with Bargo's results of same tissue conditions, it seems that there is an overall offset of 0.026 between my test results and Bargo's results in Fig.8c

Could I have a copy of source code to see how did you convert NA to acceptance angle?  such as did you assume a fiber-tissue interface during the conversion?

Owen

[Carole Hayakawa]

Hi Owen,

Thanks for posting your latest results.  Hmmm, the value 0.026 is very familiar.  I think it is specular reflection.  If you start the source in the air (in the infile in SourceInput specifying InitialTissueRegionIndex=0) and have it enter the tissue, specular reflection will occur and be subtracted.  In other words, the total light from the source = 1 = specular reflection + diffuse reflection + absorbed energy + transmitted light.  If you start the source in the first tissue layer (InitialTissueRegionIndex=1), specular reflection won't occur and total light from source = 1 = diffuse reflection + absorbed energy + transmitted light.  For air n=1.0 and tissue=1.4, the specular reflection value is 0.0277777...

Do you think this is what is accounting for the mismatch in your results compared to Bargo's?

If you'd like to see the source code, clone a copy of the code, then use something like SourceTree to check out the branch ch-141202-detector-NA-changes.  We have documentation on pulling a clone and updating to a branch for windows here:

https://github.com/VirtualPhotonics/VTS/wiki/Getting-Started-in-Windows

There are also pages for linux and os x.

Please keep me updated and let me know if you think there still might be an error in the code somewhere.

Thanks!

Carole

[future...@gmail.com]

Hi, Carole,

I may not explain the 2.6% clearly in my last post.

in my last post,  I run 2 sets of simulations, one with NA = 0.22 and the other NA = Infinite. 

Each set includes 6 tissue conditions: musp = [0.01,  0.1,  1.0,  10,  100, 1000]mm^-1, while keeping mua = 0.01 and g =0.9; n = 1.4; rho from 0 to 0.5mm;

For set 1 (NA = 0.22), total diffuse reflectance collected by rho =[-0.5, 0.5]mm, = 

[0.0001    0.0004    0.0028    0.0255    0.0462    0.0487],

For set 2 (NA = Infinite), total diffuse reflectance by rho =[-0.5, 0.5]mm, = 

[0.0004    0.0046    0.0534    0.5112    0.9013    0.9519];

then eta_c is calculated as set1./set2 = 

[15.7996    8.9972    5.3000    4.9935    5.1228    5.1191] %.

Yet, it seems that if I substract 2.6% from above eta_c,  it matches with Bargo's results in Fig.8c. at same mua;

I start the source in the air (in the infile in SourceInput specifying InitialTissueRegionIndex=0), as you can see from the attached infiles i used. 

[Carole Hayakawa]

Hi Owen,

So let me see if I understand correctly.  On your plot you level out to a nc value of 0.02419 (at x value 10^3) and in Bargo's plot 8C he lists the value as 0.0266, which is a relative difference of about 9%.  Is this the discrepancy you're pointing out?

I noticed that in your infile you specified G to be 0.9 and I think Bargo used 0.85.  Also, you specified N of the tissue as 1.4 and I think Bargo used 1.35.  So if your intent is to match his results, I would make sure the input is identical.  I'd like to point out that you can have two detectors in the same infile specifying different NA.  Just make sure that the detector Name is unique so you know which detector matches which NA.  This saves on run time.  I have attached the infile I ran that shows how to define two detectors with different NAs.

I have two suggestions:

1) make sure your input is the same as Bargo's and see if there is still a discrepancy.  If so, possibly

2) contact Bargo and compare setups.  It is possible we are missing a key element of their simulation.

With 9% relative difference with a paper I didn't write, I'm not sure I'm convinced there is a bug in the code.  Please let me know if you think there still is a bug and what we might do to debug it.

Best,

Carole

[Carole Hayakawa]

Infile with two detectors different NA.

[future...@gmail.com]

Hi, Carole,

Here is my personal thoughts.

"So let me see if I understand correctly.  On your plot you level out to a nc value of 0.02419 (at x value 10^3) and in Bargo's plot 8C he lists the value as 0.0266, which is a relative difference of about 9%.  Is this the discrepancy you're pointing out?"

---------

Not, The relative difference is not 9%, but 107.29%.   At high scattering range, Bargo's results shows that eta_c could be predicted by (NA/ntiss)^2 = (0.22/1.4)^2 = 2.47%. While my results show 5.12%, as seen from the figure in my March 26 post. 

So the percent error = ( 5.12 - 2.47)/2.47 * 100% = 107.29%.

Eta_c at high scattering range is independent from anisotropy g. At least, g-effect is negligible.  Despite other differences in simulations mechanism, the high-scattering eta_c is expected to be close to   (NA/ntiss)^2.  

Sincerely

Owen

[Carole Hayakawa]

Hi Owen,

I think I need to regroup to get a handle on what we've learned.  And just to keep it simple let's stick with the case when mus'/mua=100. 

1) On March 26, you posted a plot eta_c_1.png that shows that you determined a value for eta_c for this test case to be 5.123%.  In the associated rar file you provide, the infile is not included.  I'm assuming you used mus'=1/mm, mua=0.01/mm, n=1.4 and g=0.9.  Please verify these values.

2) On April 1, you state that Bargo's states that eta_c could be predicted by (NA/ntiss)^2.   This is a reformation of Equation (10) correct?

I'm appreciating the inconsistency you're pointing out.  Let me review the code and see if I can devise a unit test that verifies Equation (10).

Best,

Carole

[future...@gmail.com]

Hi, Carole,

"1) On March 26, you posted a plot eta_c_1.png that shows that you determined a value for eta_c for this test case to be 5.123%.  In the associated rar file you provide, the infile is not included.  I'm assuming you used mus'=1/mm, mua=0.01/mm, n=1.4 and g=0.9.  Please verify these values."

Reply: 5.123% is achieved at the case when mus' = 1000mm^-1;  mua = 0.01mm-1, g=0.9;  ntiss = 1.4;  I attached this high-scattering results in this post.

To get those 6 points in eta_c_1.png,   I run 2 sets of simulations, one with NA = 0.22 and the other NA = Infinite. 

Each set includes 6 tissue conditions: musp = [0.01,  0.1,  1.0,  10,  100, 1000]mm^-1, while keeping mua = 0.01 and g =0.9; n = 1.4; rho from 0 to 0.5mm;

set 1 (NA = 0.22), total diffuse reflectance collected within rho =[-0.5, 0.5]mm, =  [0.0001    0.0004    0.0028    0.0255    0.0462    0.0487],

set 2 (NA = Infinite), total diffuse reflectance collected within rho =[-0.5, 0.5]mm, =  [0.0004    0.0046    0.0534    0.5112    0.9013    0.9519];

then eta_c is calculated as set1./set2 =  [15.7996    8.9972    5.3000    4.9935    5.1228    5.1191] %.

---------------------------------------------------------------- 

2) On April 1, you state that Bargo's states that eta_c could be predicted by (NA/ntiss)^2.   This is a reformation of Equation (10) correct?

Reply:  Correct.

[Carole Hayakawa]

Hi Owen,

I've been working on this for a while this week and am stumped.  I found another reference of Bargo's: IEEE Journal of Selected topics in Quantum Electronics 9(2) 2003 and tried a few cases from that manuscript.  My code appears to be determining a value for eta_c that is twice what it should be.  I have reviewed the code, stepped through examples and still no luck determining where the problem is.  I think I would like to share with you the code that checks for if the photon is exiting within NA or not:

        /// <summary>

        /// method to determine whether photon direction is within NA of detector

        /// </summary>

        /// <param name="dp">photon data point</param>

        /// <param name="detectorNA">numerical aperture of detector</param>

        /// <param name="n">refractive index of region where the detector resides</param>

        /// <returns>boolean</returns>

        public static bool IsWithinNA(this PhotonDataPoint dp, double detectorNA, Direction detectorNormal, double n)

        {

            var photonDirection = dp.Direction;

            // determine if sin(theta)<=NA/n where theta is angle between photon direction and detector normal

            var cosTheta = Direction.GetDotProduct(photonDirection, detectorNormal);

            var sinTheta = Math.Sqrt(1 - cosTheta * cosTheta);

            return sinTheta <= (detectorNA/n);

            //return (detectorNA/n) >= Math.Sqrt(1 - cosTheta * cosTheta);

        }

I commented out the last line and added intermediate variables so I could view them when I'm debugging.  In the test cases we are running detectorNormal is always (0,0,-1) (positive z is into tissue).  So Direction.GetDotProduct is just picking off the Uz of photonDirection (it will have a negative Uz since escaping tissue).  I checked Scott Prahl's website because he is good about sharing code, however, the code they used for their papers is not included in his website.  As a last resort, I might contact Bargo or Scott.  I thought I'd share the code with you first to see if anything pops out at you that is wrong.  I tried several NA cases and here are the results using N=1e5 and tissue n=1.33, mus'=1/mm, mua=0.5/mm, s-d separation 2.5mm, detector width 0.6mm.

NA                           Infinity     1.33       1.0         0.39         0.22

MC eta_c                   1              1            1         0.1538     0.0476

theory (NA/n)^2         1              1            0.56    0.086       0.0274

Aside from setting NA=Infinity and 1.33, the MC eta_c values are almost twice theory.

So I think I need to enlist your help.  If nothing seems wrong with my code, I'll contact Bargo or Scott.

Thanks for any feedback you might offer.

Carole

[Carole Hayakawa]

Hi Owen,

I think I figured it out!  I was looking in the wrong place.  The detector NA code is fine.  The way the code works is that we process the photon until it escapes the surface and have it refract into air.  We do this for cases in which the photon is captured by a camera above the tissue for instance.  If you notice in my prior post, all of the "MC eta_c" values look goodI if the theory uses n=1.  So setting FinalTissueRegionIndex=0 works fine.  However, I read the Bargo papers more closely and they assume the detector is sitting right above the tissue so I am assuming the photon does not refract before it is detected.  This would be the case when FinalTissueRegionIndex=1.  The code needs to determine the direction prior to refraction into air and this code is in question.  I did a quick fix to the code to run this case and the results look good and agree with Bargo.

It will take me a little bit to fix the code correctly and add unit tests to check for this change to the code.

I'll keep you posted.

Best,

Carole

[future...@gmail.com]

Hi, Carole,

Glad to hear it.

I recommend you set a higher value for mus' in your simulations, such as mus' = 100mm^-1;   since eta_c = (NA/n)^2 is only for high -scattering case. 

Best

Owen

[Carole Hayakawa]

Dear Owen,

Upon working on a correct code (no quick fix), I realized that Bargo et al. simulation is designed differently then ours.  I check with Scott Prahl and he is fairly certain that the code that they used applies a Fresnel check at the tissue/detector interface differently than at the tissue/air interface.  In other words, when photons arrive a the tissue surface their code checks Fresnel using n=1.33 (tissue) vs n=1.45 (detector) when under the detector and reflects/refracts accordingly and checks Fresnel using n=1.33 (tissue) vs n=1.0 (air) when not under the detector.  Our code treats the tissue/air interface the same across the entire surface.  So using our code as it currently stands will never produce their results. 

The quick fix that I tried is that I specified the detector NA to be (NA/n_tissue) equivalent.  So in their paper, if they specified NA=0.22 for a tissue with n=1.33, then the "effective" NA in air would be (0.22/1.33)=0.1644.  Using this value for the detectors NA in the infile produces results that are closer to theirs.

I'm going to discuss the possibility of modifying the code to mimick Bargo's, however, this change would couple the detector definition with the photon transport, something that is currently not supported in the code. 

So as it stands now, I feel our code is correct for detector definitions within the air

Please let me know any questions or comments you might have.

Best,

Carole

[future...@gmail.com]

Hi, Carole,

I find a potential issue in your code posted in April 4, which may explain why your quick fix works. "So in their paper, if they specified NA=0.22 for a tissue with n=1.33, then the "effective" NA in air would be (0.22/1.33)=0.1644.  Using this value for the detectors NA in the infile produces results that are closer to theirs."

In your code,  return sinTheta <= (detectorNA/n);

 "n" is refractive index of region where the detector resides</param>,   Since detector resides in air,  so n = 1.0

However, we calculate the sine of acceptance angle in tissue  by  sin(AcceptanceAngle) = (detectorNA/ntiss),  we need use refraction index of the region which the detector fiber directly contact.

In other words, The code itself consider n = 1.0 when calculating acceptance angle. So, this explains why when user specify an effective NA = (NA/ntiss) could produce a result closer to Bargo's.

You can trace your code to see how is n defined.

I tried one case [musp = 100mm^-1, mua = 0.01,  g = 0.9, ntiss = 1.4, NA = 0.22], the resulted collection efficiency = 2.58%, compared to 2.47% from Bargo. Percent error 4.5%.

------------------------------------------------------------------------------------------------------------------------------------------

 /// <summary>

        /// method to determine whether photon direction is within NA of detector

        /// </summary>

        /// <param name="dp">photon data point</param>

        /// <param name="detectorNA">numerical aperture of detector</param>

        /// <param name="n">refractive index of region where the detector resides</param>

        /// <returns>boolean</returns>

        public static bool IsWithinNA(this PhotonDataPoint dp, double detectorNA, Direction detectorNormal, double n)

        {

            var photonDirection = dp.Direction;

            // determine if sin(theta)<=NA/n where theta is angle between photon direction and detector normal

            var cosTheta = Direction.GetDotProduct(photonDirection, detectorNormal);

            var sinTheta = Math.Sqrt(1 - cosTheta * cosTheta);

            return sinTheta <= (detectorNA/n);

            //return (detectorNA/n) >= Math.Sqrt(1 - cosTheta * cosTheta);

        }

------------------------------------------------------------------------------------------------------------------------------------------

[future...@gmail.com]

Carole,

I attached the raw simulation result in this post. I FinalTissueRegionIndex= 1 in the infile.

[Carole Hayakawa]

Dear Owen,

Thank you for your posts and data.  The code that you are using processes the photon until it escapes the air, that is, the final direction is the refracted direction into air.  So using that direction only makes sense with checks using n of air.  For infiles that specify FinalTissueRegionIndex=1, I "rewind" the photon to the direction prior to refracting into air and use the n of the tissue.  However, these two checks are the same since n1 sin(theta1) = n2 sin(theta2).

I *do* appreciate your input and have discussed with my colleagues the need to have code that allows the user to define a detector with a refractive index n and NA *adjacent* to the tissue surface.  This would entire changing the transport of photons that hit the surface under the detector and if they refract, tally them if they refract into the NA and then kill them.  Since that is a situation that occurs in experiments often, we feel it would be good to have coded.  It would also allow you to run simulations equivalent to Bargo's.

I plan to work on this code on this branch you and I have been sharing.

I'll keep you posted.

Carole

[future...@gmail.com]

Hi, Carole,

Please do not hesitate to contact me if you need me to do some test on your new code.  I am glad to help.

Currently I am digesting your code on the branch you have been sharing with me, hoping to understand it better to provide better help.  I has programming experience but new to C#. Any illustration on your code and any instructions on how to get on track quickly will be greatly appreciated.

Best.

Owen.

[Carole Hayakawa]

Hi Owen,

Thank you very much for your great offer to help.  Once I have a plan in place, I will do the work on the branch.  I'll definitely update to this discussion when I do.

Unfortunately we don't have instructions about how to come up to speed with the source code.  If you have the source code cloned locally, the Vts.Desktop is the library that contains all of the MonteCarlo classes.  The Vts.MonteCarlo.CommandLineApplication is the program and with the Vts.Desktop is what comprises the MCCL download.

I will post again soon.

Best,

Carole

[Carole Hayakawa]

Hi Owen,

In the Issues tab I've been posting my code updates for this modification.  I've created a new zip file and added in a sample infile named infile_surface_fiber_detector.txt.  I ran this infile with 1e5 photons and compared the results of NA=0.22 vs NA fully open and getting an eta_c value of 0.0291, however, if you take into account the standard deviation of the two estimates, it encompasses the value 0.0266 that Bargo generates.  

I gave this an unofficial version of 0.0.1 so that it is different from the prior zip file.  The link is:

https://drive.google.com/file/d/1U-M47c3LhW-PHqJ3jFeGLKBiPGEIbTvG/view?usp=sharing

Please let me know how the testing goes and if you think this isn't an improvement over using the prior code and an effective NA.

Thanks for your help testing!

Best,

Carole

[future...@gmail.com]

Hi, Carole,

Excited to hear it!

But the link you shared still points to previous version MC_v0.0.0Beta.zip that was built on 3/8/2019. 

[Carole Hayakawa]

Hi,

Let me try again.  How about this?

https://drive.google.com/file/d/13jAzjOEvGijE3K_eLgMoRSPz8EXO62f7/view?usp=sharing

Best,

Carole

[future...@gmail.com]

Hi, Carole,

Please check the attached RAR file for my test results, figures, and raw data as well as infiles. 

At high scattering condition mus' = 100mm^-1, the collection efficiency = 5.2% with stdev = 0.02%.   Bargo's formular = (NA/ntiss)^2 = (0.22/1.4)^2 = 2.47%. 

To my understanding about his paper, his formular above is only valid for high scattering condition.  The case mus' = 1mm^-1 you tested, may not be considered as high scattering condition.

Owen.

[Carole Hayakawa]

Hi Owen,

Thank you very much for these results.  After reviewing your code, I have a few comments:

1) The first detector that you set up in all cases sets the detector N to be 1.0, that of air, however the FinalTissueRegionIndex is set to 1 (tissue) so not sure what to make of what the results should be. The code for this detector is new and does not conform to the code in the R(rho) detector that we were analyzing previously.  So to compare with Bargo, I would not consider this case.  If you'd like to see how our prior results compare, include the ROfRho detector with the specifications were setting above.

2) The eta_c2 in your code is the calculation I would focus on (blue line in your plot).  Looking at Bargo Figure 8(c) and his plot for mua=0.1/cm and mus'=10/cm (which would agree with your infile mua=0.01/mm and mus'=1/mm), your calculated value of eta_c=0.0273 looks like it might agree with his.

3) You then ran a test case with mua=0.01/mm and mus'=10/mm.  Here the blue plot starts to rise and in Bargo's plot, the plot looks like it is monotonically decreasing to reach the asymptote 0.0266.  Your reflectance open Ropen = 0.2199 +/- 0.0001 (1 SD or 1 sigma) and for NA=0.22 R0p22 = 0.0086 +/- 0.00005.  So even accounting for a possible 3 sigma error in the mean value of both reflectances, we still cannot obtain a value of eta_c close to 0.0266.

I ran a test case of my own for mus'/mua=1000 and get a eta_c value of 0.036.  Your value is 0.0392.  So I understand your question about why this is not agreeing with Bargo's results.  It is curious that our eta_c values look good for mus'/mua=100 but not mus'/mua=1000.  I can't think of anything obvious in the code that would cause this.  Looking at your results for mus'=1 NA0p22.txt shows a TallyCount of 2957 and Open shows a TallyCount of 107863 producing a ratio of 0.027.  For mus'=10 NA0p22.txt shows a TallyCount 27611 and for Open TallyCount 700100 giving a ratio of 0.039.  It is kind of surprising that these ratios match the eta_c values so closely.  But that gives me a hint about where to look for a possible bug in the code.  I will proceed with that.

Please let me know if you disagree with any of my comment above.

Thanks again for your help in debugging the code. 

Best,  Carole

[future...@gmail.com]

Hi, Carole,

1). I agree.

2). I agree.

3).  In my infiles, I specified refraction index N = 1.4 for tissue, not 1.35 as Bargo used. That is why at high-scattering region, eta_c = (0.22/1.4)^2 = 2.47% for my simulated data, according to bargo's formula.

4).  I agree.  

[Carole Hayakawa]

Hi Owen,

I thought I would just give you status.  I have been working on debugging the code *steadily* since my last post.  I think I ran about 30+ debug cases to try to figure out the bug in the code.  I found very minimal code changes that improve results for a more general input, but won't change your cases much.  I then returned to the theory in the paper.

I stared at Eqn. (10) and realized the "cos(theta)" term might mean that possibly I'm suppose to multiply the photon weight by its Uz prior to adding it to the reflectance mean.  This kind of made sense since for a fiber, the light projected onto the normal might be needed. I wrote a little code that selects Uz ~ unif(0,1) and then calls my method to determine whether it is within NA.  If it is I add Uz to a sum for NA=0.22 and NA=Inf independently. I ran 1e6 random Uz and when I take the ratio I get 0.0247.  It also produces the correct eta_c value for NA=0.39.  So this told me that when I tally (if in NA), I should multiply photon weight by its Uz.  I coded this, however it worsened my results.

But I'm feeling Eqn. (10) now.  It makes sense that as mus' is increased, the eta_c values converge because as mus' is increased, l* is shortened and the light field around the source becomes isotropic faster spatially.  I am missing something else in my code.

I'll keep at it.  I reached out to Scott Prahl again yesterday to see if he any suggestions for me. 

I'll keep you posted.

Best,

Carole

[Carole Hayakawa]

Hi Owen,

I *think* I have it.  I returned to some archived C code and was able to quickly code up what I thought the code should be doing.  It gave me results that look promising.  So I returned to the C# code and coded up something equivalent.  The new design is a bit cumbersome because in the infile you have to specify a TissueRegion representing the fiber on the surface that requires a Center, Radius, and OpticalProperties which will specify N.  These characteristics need to match those specified in the SurfaceFiberDetector.  I plan to write code to validate input so that it catches if those are inconsistent, but I wanted to get the code to you so that you can see how it works for you.

I ran N=1e5 photons, with mus'=10/mm, mua=0.01/mm, n=1.4 of tissue, n=1.4 in fiber, NA=0.22 and obtain an eta_c=0.0262 and using mus'=100/mm I get eta_c=0.0277.  I would be more happy if the second value was closer to 0.0247, but possibly I need to run more photons.

I uploaded a new zip file with version 0.0.2 in google drive with link:

https://drive.google.com/file/d/1bS6ahXGB1tayKyUoHTW7UQRGI1drbu8z/view?usp=sharing

I hope this improves the results for you.  Please let me know if you feel there is still something wrong with the software.

Thanks for your help!

Best,

Carole

[future...@gmail.com]

Hi, Carole,

Glad to hear! I will test it later and post here.

[Carole Hayakawa]

Hi Owen,

I hope you are well.  I was just wondering how this interesting study you were pursuing ended up.  Did you come to any resolutions concerning the results described in the literature and those obtained by the MCCL?  I'd be very interested to find out if you feel there are still discrepancies.

Thanks,

Carole

[Owen Sun]

Hi Carole,

Sorry for my looooong delayed reply. I did not get time to test V0.0.2, and I adopted an empirical formula of collection efficiency from the following paper for my project. 

Please go to page 6/16, and check equation (4). It may be useful for comparison with results from V 0.0.2 

Kanick, S., et al. (2011). "Measurement of the reduced scattering coefficient of turbid media using single fiber reflectance spectroscopy: fiber diameter and phase function dependence." Biomedical Optics Express 2(6): 1687-1702.

The above formula is empirically-obtained in steady-state monte-Carlo simulation. In other situations, e.g. the time domain, the results may be different. 

I believe it is of great value for MCCL to inherently have an NA feature in its detector. Hope to see feature in future releases.

best,

Owen

[Carole Hayakawa]

Dear Owen,

Thank you for your reply!  I will check out the paper you reference.

We now have detector NA in all of our detector definitions.   This code used to be on a branch.   It is now part of the master branch and included in recent releases.

Best regards,

Carole