This project has moved. For the latest updates, please go here.

Using command-line MC for obliquely incident Gaussian beam source with variable angle and beam waist in layered media

Sep 14, 2012 at 3:04 PM

Can the command-line MC be used to model reflectance as a function of (x,y) from layered tissue due to a Gaussian beam with a variable beam waist and incident from a variable angle? Can you help me define these conditions within your current format?
Also, I know when we were visiting last year we used Mono on my mac to run the command-line in Terminal (and I believe we got it working then), should that still work with the current version?
Thank you!

Coordinator
Sep 14, 2012 at 5:12 PM
Edited Mar 13, 2013 at 5:28 PM
Hi Shelley,

I know Carole is working with you on the Monte Carlo, I just wanted to respond about running the command-line application using Mono on the Mac. Yes It should still work with the current version, let me know if you have any problems though and I will troubleshoot.

Thanks,
Lisa
Developer
Sep 14, 2012 at 7:29 PM

Hello Shelley!  Great to hear from you.

We have a ROfXAndYDetectorInput that specifies collecting reflectance from a x-y grid on the surface of the tissue, R(x,y).  The two inputs are the x range and y range.  In the Monte Carlo Command Line Application download there is an example of this detector input infile infile_one_layer_all_detectors.xml.

Now for the source.  There is an infile in the download named infile_Gaussian_source_one_layer_ROfRho.xml that specifies a Gaussian source.  I think you might have used this infile last time?  Anyway, the source is called CustomCircularSourceInput.  If you go to the write-up on

http://www.bli.uci.edu/vp/wiki/index.php/Conventional_Monte_Carlo

under "Sources:....click *here*", there is a description of the Gaussian source starting on page 2.

To angle the Gaussian, there are two ways to angle the source.  Look on pages 4-5 of the write-up.  The first way is item 2 on page 4 in the write-up.  This samples from the Gaussian beam that would be normal to the tissue BUT then rotates the photon's direction vector by (theta,phi).  This is performed by setting the 8th parameter to CustomCircularSourceInput and specifying the (theta,phi) pair that you would like to rotate to.  The second way is item 3 in the write-up.  This rotates the entire source towards a given direction set by the 6th parameter to CustomCircularSourceInput.  This changes the initial starting position of the photons.

Let me know if you have questions on any of this or if you have trouble piecing together an input file with your desired source, tissue and detector definitions.

Thanks for your continued interest in our software!

Carole

Sep 17, 2012 at 10:46 PM

Thank you both! I think I got it running on my Mac over the weekend. Your descriptions are quite helpful! :)

Feb 28, 2013 at 9:03 PM
Hello! I have a, perhaps odd, additional question: When I run MC_V1.0.5Beta with infile_Gaussian_source_one_layer_ROfRho.xml (parameters modified for our problem) - I'm curious why the reflectance reduces to 1 away from the origin instead of 0 - any ideas? I haven't noticed this behavior in the MC on the website VTS. Thank you!
Developer
Feb 28, 2013 at 9:43 PM
Hi Shelley. I just ran the infile_Gaussian_source_one_layer_ROfRho.xml for the latest download MC_V1.0.10Beta and plotted the R(rho) results (for 100 photons). It appears that the reflectance is going to 0 as rho increases. I'm using the default infile parameters though. Note that the last bin of R(rho) stores all photon weight that escaped the surface beyond the last bin, so that bin might show a rise in the reflectance from the second to the last bin. We do this so that if you sum over all rho bins it will equal the total diffuse reflectance. Can you send me your infile and possibly a plot of your results? Do you think it is particular to V1.0.5Beta or do you get the same results for the latest download?
Carole
Mar 1, 2013 at 12:39 AM
I'll install the new version and see if that fixes it (mine is kinda old). It's not a bump, it looks like the reflectance is converging to 1 as rho increases - which seemed odd. If the new version doesn't fix it, I'll send you my files & plot - Thank you!
-Shelley

On Thu, Feb 28, 2013 at 2:43 PM, hayakawa <notifications@codeplex.com> wrote:

From: hayakawa

Hi Shelley. I just ran the infile_Gaussian_source_one_layer_ROfRho.xml for the latest download MC_V1.0.10Beta and plotted the R(rho) results (for 100 photons). It appears that the reflectance is going to 0 as rho increases. I'm using the default infile parameters though. Note that the last bin of R(rho) stores all photon weight that escaped the surface beyond the last bin, so that bin might show a rise in the reflectance from the second to the last bin. We do this so that if you sum over all rho bins it will equal the total diffuse reflectance. Can you send me your infile and possibly a plot of your results? Do you think it is particular to V1.0.5Beta or do you get the same results for the latest download?
Carole

Read the full discussion online.

To add a post to this discussion, reply to this email (virtualphotonics@discussions.codeplex.com)

To start a new discussion for this project, email virtualphotonics@discussions.codeplex.com

You are receiving this email because you subscribed to this discussion on CodePlex. You can unsubscribe on CodePlex.com.

Please note: Images and attachments will be removed from emails. Any posts to this discussion will also be available online at CodePlex.com


Jun 12, 2013 at 8:08 PM
Hello,

I've been trying to debug the issue with MC on my Mac, to no avail. I updated from my previous version of matlab (2007a) to the current one (2013a) and now seem to be getting a whole new problem with the reflectance plots. I'm going to send a .zip with relevant documents (plot, load_results_script.m, the folder from the MC run, and the infile). I can't seem to attach it here.
I ran it with the latest MC download version on my Mac, and the only edits I made to the downloaded docs are the #of photons in the infile, and the path, "datanames", and "outdir" in the .m file.

Any ideas about what's happening?

Thank you!
-Shelley
Developer
Jun 14, 2013 at 9:21 PM
Hi Shelley,

I've executed the infile you sent me on linux with download version MC_v1.0.10Beta and my results and yours look identical. Here are yours. Here are mine. Now as I mentioned above, the reason the plot rises at the last rho is because we tally all photons that escape the tissue surface beyond that bin, in that last bin. Note that the Gaussian beam source default OuterRadius is set to 3mm as it is in your infile so that presumably gives the reflectance plot a bit of a "hump" near rho=0.

What problem with the reflectance plot were you concerned with?

Carole
Jun 14, 2013 at 10:26 PM
Hi Carole,

Thank you! I was concerned because it is so negative, and when I ran previous infiles I was experiencing results that appeared to be too positive (here it appears to be approaching -4, and in my previous runs it was approaching 1). When I saw the negative result, I thought it there was still an issue because I don't usually see that happen with our reflectance model unless there is an issue.

The spike at the end does not concern me, but the overall values of the reflectance were my concern.

So, I shouldn't be concerned?

Thank you!
-Shelley


Developer
Jun 14, 2013 at 11:12 PM
Hi Shelley,
I don't have much experience with results from a Gaussian beam. But if you use the same optical properties (mua=0.01/mm mus'=1/mm g=0.8 n=1.4) and run the Forward Solver Engine=NURBS using our GUI and even though this is a point source, you'll see that reflectance y-axis range is similar (from -1.5 to <-3.5 on a log scale). So I would believe your results. You could also squeeze the size of the OuterRadius of the source and see how that affects the drop off rate near rho=0.
Anyway, I don't think you need to be concerned. Remember though that if you start angling the source, then you'll need to switch to using ROfXAndYDetector.
Let me know if you run into any other problems.
Carole
Jun 15, 2013 at 1:18 AM
Thank you Carole! I appreciate your help!


Developer
Mar 21, 2014 at 8:21 PM
I am posting an email exchange SRohde and I had that continues this topic:

SRohde asked: For an obliquely incident Gaussian source, is there a limitation on the range of angles for which it is accurate?

hayakawa response: Whatever angle you choose, the results should be valid, HOWEVER, the variance in the results may vary greatly.  So you could turn on "tallySecondMoment" in SimulationOptions and use our matlab post-processing code to determine the standard deviation in your results.

janakarana added: I would like to add something to Carole's reply.
Dark horizontal line after the rotation shows your actual tissue layer boundary. We assume original source is normal and it points downward. In item 2, the tissue layer boundary is horizontal. Then, the source is rotated by a given angle pair. Let me give some examples and their angle pair settings for your reference.
        1. rotate pi/6 counter-clockwise around x-z plane  ----> set angle pair to (pi/6, 0). 
        2. rotate pi/6 clockwise around x-z plane ----> set angle pair to (pi/6, pi). 
        3. rotate pi/6 counter-clockwise around y-z plane  ----> set angle pair to (pi/6, pi/2). 
        4. rotate pi/6 clockwise around y-z plane ----> set angle pair to (pi/6, 3pi/2). 
In item 3, you can see that tissue layer boundary is not horizontal. If someone has a tilted tissue layer and wants to use a normal incident, he/she has to use this rotation. This feature is included for our future developments and it is not applicable for horizontal tissue layer boundary.

Since your tissue layer boundary is horizontal, you can use the setting stated in item 2.  The angle limit for theta is pi/2. If theta >= pi/2, photons will not go into the tissue. As Carole mentioned, the results are valid for any theta (< pi/2).

SRohde asked: I'm still having issues with an angle just over pi/4, I'm actually receiving zero reflectance for this case. I extended the values of x and y in an attempt to ensure they would capture reflectance, but am receiving zero - any ideas why this might be happening?

janakarana response: This is strange. I have to check the code. Could you give me your source settings that provides zero reflectance. What is your beam radius and how many photons do you simulate?

SRohde response: I'm using 1000000 photons generally, but have tried 10000000. With a beam radius of FWHM = 5, outer radius then set to 15. Theta = 0.815868214585956

Let me know if there are any other parameters you need. Here is the output.xml file from the run, if that is helpful. (I did not use the Fluence for anything, only ROfXAndY)

janakarana response: Shelley mentioned that she was receiving zero reflectance for oblique angles over 45degree. I checked the code for customcircularsource to see whether any problems with the calculation of direction and position. For the given theta (0.8158...), direction coefficients are correctly calculated; ux =0.72832, uy=0, and uz=0.685236. The position values lie within the range. The direction and position values are properly passed to photonDataPoint parameter "DP". Then Snell's law is applied to calculate the correct direction. I don't see any problem up to that point.

It is strange to see zero reflectance. Shelley included her infile in her email. Would you able to do a test run and see? I am copying Shelley on this email.

hayakawa response: I found the bug.  ROfXAndYDetector was not incrementing TallyCount.
  The fix is to add the line:
  TallyCount++;
  after the code section:
  if (_tallySecondMoment)
  {
    SecondMoment[ix,iy]+=photon.DP.weight * photon.DP.weight;
  }

  so it looks like:
  if (_tallySecondMoment)
  {
    SecondMoment[ix,iy]+=photon.DP.weight * photon.DP.weight;
  }
  TallyCount++;

  Let me know if this works for you or not.  In the meantime, I'm going to determine why my unit tests did not find this bug :(

  I'm sorry for your problems with our code Shelley.
  Thanks for helping us improve our code base.
  I will update the software, add unit tests validate the fix and this fix will be in the next MC release.
Developer
Mar 21, 2014 at 8:40 PM
I just realized that the bug I found is a true bug but may not the root of your problem. I removed the bug fix and I am still able to plot a non-zero R(x,y) using the matlab post-processing software in matlab/post_processing/monte_carlo/simulation_result_loading/.

So let me know if you're still having problems or have any questions.
Mar 21, 2014 at 10:06 PM
Hi Carole,

I knew there was a problem with those runs because they would take a couple minutes for 1000000 photons, whereas I'm used to those runs taking close to a day for other angles.

So, you don't think that line will fix the problem?
I didn't think the zero reflectance was tied to the post-processing because of the time each run was taking.

You know the code better than I do - please help me understand.

Thank you!
-Shelley


Developer
Mar 21, 2014 at 10:49 PM
Try this: Use the infile that you previously sent Janaka and make sure you're using the Monte Carlo Command Line Download 1.0.12. Change the number of photons to 1000. You have the seed set to 0 in simulation options, so you and I should reproduce the same results. Now run this infile (it should take about 6 seconds). A folder named oblique_Gaussian_two_layer_kienledata_ep16_rot is created. In the MC download there is a file called load_results_script.m. Open that file and edit line 11 to be "datanames = { 'oblique_Gaussian_two_layer_kienledata_ep16_rot' }; Now run this matlab script. You should see a R(x,y) plot that has 931 reflected photons (I know this amount from running the code with the "TallyCount++" addition). Let me know how it goes. Carole
Mar 21, 2014 at 11:08 PM
I still have 0 reflectance, and the MC runs in under 1/20th of a second :/ There must be something wrong on my end. This doesn't make sense.

I reinstalled v1.0.12, used its mc.exe, and it still gave me the same problem. I'm not sure how that's possible when we're using the exact same .xml file... unless it's because I have to run tests on my Mac right now - but I had exactly the same problem on the Windows box I use at home. I'm completely puzzled by this. Any ideas?

Thank you for all your help!


Developer
Mar 21, 2014 at 11:25 PM
How about this: Run the infile_one_layer_all_detectors.xml file that is included in the download (takes about 3 seconds). After running, substitute sub-folder name "one_layer_all_detectors" on line 11 of script load_results_script.m and run in matlab. You should see a series of plots with non-zero results, including Figure 7 which is R(x,y).
Mar 22, 2014 at 8:02 AM
here's the output from the processing after that run (none of which are explicitly zero the way my infile is):
Total reflectance captured by RDiffuse detector: 0.5841
Total reflectance captured by ROfAngle detector: 0.27176
Total reflectance captured by ROfRho detector: 0.59366
Total reflectance captured by ROfRhoAndAngle detector: 0.27967
Total reflectance captured by ROfRhoAndTime detector: 5.9366
Total reflectance captured by ROfXAndY detector: 0.005841
Total reflectance captured by ROfRhoAndOmega detector: 0.59366
Total transmittance captured by TDiffuse detector: 8.4416e-17
Total transmittance captured by TOfAngle detector: 1.7537e-16
Total transmittance captured by TOfRho detector: 1.3503e-17
Total transmittance captured by TOfRhoAndAngle detector: 2.8051e-17
Total absorption captured by ATotal detector: 0.3859
Absorbed energy captured by AOfRhoAndZ detector: 4.2162
Fluence captured by FluenceOfRhoAndZ detector: 421.6182
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfRhoAndZAndTime detector: 4207.7669
Fluence captured by FluenceOfXAndYAndZ detector: 192.9478
Radiance captured by RadianceOfRhoAndZAndAngle detector: 60.4136
Radiance captured by RadianceOfXAndYAndZAndThetaAndPhi detector: 12097.2535

However, none of the other (smaller) angles I used seemed to have any problem. For some reason, it's just that one angle causing me grief. The fact that it's less than pi/2 gave me concern about the rest of my results - but the other results match our model well, so I had no clear reason to worry aside from the angle I sent you.