Mar 13, 2012 at 9:57 AM
Edited Mar 13, 2012 at 9:59 AM

Hi:
First of all thanks so much for this amazing open source software. I have a quick question. I am currently
trying to simulate the light transmission and reflectance in a turbid medium , so I am using the Monte Carlo solver.
At the end of calculations, when I am trying to open the map view, nothing appears. Is that because map view
option is not implemented yet in MC solver? Also can I get the phase function versus angle?
Thanks for your time.
Alireza




Hi Alireza!
Thanks for the question. If you're using the GUI on the web, you're right  visualizing internal fluence for Monte Carlo simulations is not yet implemented in the online tool.
If you save the results to your local machine (click "Save Simulation Results"), you can read the doubleprecision binaries with the tool of your choice. For example, you can use the free ImageJ tool using Import>Raw>64bit Real. If you have Matlab,
we actually already supply some simple scripts for plotting and visualizing the results to get you started. You can find those tools in
this source folder, or as part of the desktop simulator download
here.
Please let us know if this works out for you, or if we can offer more help!
Yours,
David




Forgot to answer your second question. For now, we only allow two types of phase functions: HenyeyGreenstein (described in
this article) or Bidirectional, enumerated
here. The default is HenyeyGreenstein (Bidirectional is only really useful for performing testing, i.e. comparisons with analytic solutions).




Dear Alireza,
Thank you very much for your interest in our project!
When using the Monte Carlo Solver Panel, you need to specify an input file that generates map type data, for example, fluence.
In our download of prototype input files, there is one called "infile_one_layer_ROfRho_FluenceOfRhoAndZ.xml". If you load this input file and run, there should be fluence results in the "Map View" and spatiallyresolved reflectance results in the "Plot
View".
The phase function in the Monte Carlo simulation specified in the prototype input files is the HenyeyGreenstein phase function. In these input files, an anisotropy coefficient of g=0.8 is specified, but can be modified for your particular problem.
We currently don't have the capability to plot the phase function versus angle, however, it can be determined using the formula described on [http://en.wikipedia.org/wiki/Monte_Carlo_method_for_photon_transport]. Note that we also have the capability
to specify (via the input file) bidirectional scattering and in the future plan to add Mie phase functions.
Please let us know if you still have any problems or questions.
Thanks,
Carole




Sorry, disregard my first comment re: viewing fluence. Carole's right, we do have Map View enabled for the specific case FluenceOfRhoAndZ detector, as she mentions above. David




Hi:
Thanks really so much for all the comments. Very helpful and I could see the fluence versus Z. But is there any possibility to see the reflection as well (Z<0)
So basically not onlly visualizing the internal but external scattering as well. And at the end what I want to do is looking at the total reflection, which is the integration of reflection over back hemisphere, Can you suggest me a efficient way
to do that?
Thanks again in advance.




Never mind, everything is clear. Just a quick question, I also tried to take a look at results in imagej but I could not load the files.
Just to make sure about the steps: I should save simulation results, then importing the zip file in raw format (64bit real) !
Thanks for your help.




Hi Alireza,
Sorry, I was a little ambiguous. Here are the full steps I used to do this in the Monte Carlo Solver Panel:
1) Click "Download Prototype Input Files". Save zipped file with a name such as "inputs.zip"
2) Unzip "inputs.zip" to a local "inputs" folder.
3) Click "Load Input File". Select "inputs/infile_one_layer_ROfRho_FluenceOfRhoAndZ.xml".
4) Change Number of Photons to 1000
5) Click "Run Simulation" and wait for simulation to complete
6) Click "Save Simulation Results". Save zipped file with a name such as "results.zip"
7) Unzip "results.zip" to a local "results" folder
8) Open ImageJ and use "File>Import Raw..." to select "results/FluenceOfRhoAndZ" (the binary, not the matching xml file)/
Import>Raw Settings are:
Image type: 64bit Real
Width: 100
Height: 100
Offset to first image: 0
Number of images: 1
Gap between images: 0
White is zero: false
Littleendian byte order: true
Open all files in folder: false
Use virtual stack: false
9) Chose "Process>Math>Log"
10) Chose "Image>Adjust>Brightness/Contrast" (I use Ctrl+Shift+C), and click "Auto"
You should see this: http://i.imgur.com/9gsq4.png



Mar 14, 2012 at 8:14 PM
Edited Mar 14, 2012 at 8:14 PM

Just to follow up, here are the matlab commands to do the same thing:
0) Run steps 17 above
8) Navigate Matlab to the "results" folder
9) Run the following:
h = fopen('FluenceOfRhoAndZ','rb');
fluence = fread(h,[101 101],'double')';
fclose(h);
figure;
imagesc(log(fluence));
colormap gray
Note that ImageJ and Matlab report NaN's as white and black, respectively (for cases where we took the log of 0...voxels deep beneath the source that had no weight with a 1000 photon run)




Hi:
Thanks really so much. I could go through all the steps and find what i was looking for.
Can you please refere me to a place to find the exact definition of fluence and R. Because
at the end I want to calculate the amount of light which is backscattered and probably I need
to do some integration and post processing.
thanks so much for all your help.
Alireza




No problem. We actually have a detector for the whole diffuse reflectance (no spatial binning)  it's called RDiffuse and it's the first one listed in the template file "infile_one_layer_all_detectors.xml"
Our documentation is still catching up  you can read a little about some of the detectors
here. Generally our units are in mm, ns, and radians. Some of the detectors you might be interested in are:
RDiffuse [unitless fraction]  fraction of photons hitting the tissue that enter the tissue and diffusely reflect. Should be <1 and equal to (1RSpecular)
TDiffuse [unitless fraction]  fraction of photons hitting the tissue that enter the tissue and transmit.
ROfRho [1/mm^2]  Spatiallyresolved diffuse reflectance in radial bins. Units are a fraction of the light input per unit area detected (e.g. the surface area of the detector)
ROfRhoAndTime [1/(mm^2*ns)]  Same as ROfRho, but also timeresolved, thus the extra "pertime" units.




Thanks again. two questions: when you are saying fraction of photons hitting the tissue that enter the tissue and diffusely reflect, are you accounting for fresnel reflection from the surface as well as a refractive index miss match and just to make sure,
i am assuming Integral of ROfRho [1/mm^2] over the entrance area is RDiffuse.
Thanks,
Alireza




Yes, we account for Fresnel and index mismatch. We actually sample the Fresnel reflection stochastically, just like the rest of our processes.
In the default MultiLayerTissue, we also handle index mismatch between layers. Internal mismatch is not currently implemented for SingleInclusionTissue, however.
And yes, ROfRho should integrate (using an axisymmetric Riemann sum) to RDiffuse. We include photons out of the rhorange of the ROfRho detector in the final bin, so as to capture all the data. In the future, we hope to have this behavior as an option in
the input file.




Hi:
So I could be able to do the simulations. But I have a concern:
I am modelling a single layer (airtissueair), and I am using g=0.6, mus'=20/mm, n=1.5.
What I want to capture is the total amount of light get reflected. I was changing the thickness of the tissue from 0.1mm to 10mm.
I am not observing a change in ROfRho. Does it make sense from physical point of view?
Thanks,
Alireza




I forgot to mention that the mua=0.




Hi Alireza,
Sorry you're having trouble. I just ran the simulation above through our online GUI using the prototype infile: "infile_one_layer_ROfRho_FluenceOfRhoAndZ.xml".
I ran 10,000 photons and modified the singlelayer thickness from 10mm (blue) to 0.1mm green. I definitely see differences (screenshot
here) between the two that make general physical sense (lower reflectance for the thinner layer, especially at longer sd separations)
. In the screenshot, I've plotted in "loglog" form to highlight the differences close to the source.
Can you reproduce this plot? If so, can you send us a repro to get the equivalent results you've calculated?
Thanks,
David




Hi David:
First thanks so much for the response. I did not expect an answer in weekend and sorry for asking questions.
I exactly did the same thing and got the same curves. http://i.imgur.com/oDzTi.jpg
The problem I am trying to solve is: Abeam of laser is passing through a tissue with thickness of 0.1 to 10 mm.
I want to observe the reflection, which I could, and now I am trying to look at the radiance in different angles (pi/2, pi/2)when I am
looking at the tissue from below. So doing that, I used the command line and used all detectors module and looked at the RadianceOfRhoAndZAndAngle.
Is it the right way? Also I do not understand the plot as well. Which axis is the angle? Where is the place that we measure angle from? and why is it a function of rho and Z?




Hi Alireza,
It sounds like you are interested in the angleresolved transmittance and the detector TOfRhoAndAngle. The angle here is the polar angle measured from the positive z axis. Positive z is perpendicular to the tissue surface, in downward direction. So
this detector checks the direction cosine in the zdirection, uz, of the transmitted photon and takes the acos of that to bin. For transmittance, 0<=uz<=1 and so 0<=cos(uz)<=pi/2. Keep this in mind when you specify your angle
bins for the detector.
RadianceOfRhoAndZAndAngle is meant to be a detector *internal* to the slab where z specifies the depth of the surface measuring radiance. This detector detects photons traveling through this surface in a downward direction (uz>0). Again the
angle is the polar angle.
Carole




Hi:
Thanks for the reply. I tried to do that and I did not find a place to set the angle bins for the detector. Here is a screen shot of
my TofRhoandAngle detector with parameter. http://i.imgur.com/9jH8L.jpg
Also I attached the Matlab results. Could you please let me know where is the problem. And also I want to know where we are placing the detector, and why is its coordinate is only a function of Rho and not Rho and Z for example? Btw, I am assuming Rho is
Rho=sqrt(x^2+y^2) with z normal to tissue surface.




From your matlab picture, it looks like you've only got *one* angle bin specified (the default value in infile_one_layer_all_detectors.xml). But from your XML Notepad screen shot it looks like you specified 101 angle bins ('Count' under the 'Angle'),
so I'm not sure about whats happening there. Note that TOfRhoAndAngle is a transmittance detector so it detects all photons exiting the bottom of the slab. The detector has circular rho bins (so your assumption on rho is correct). The Z of
detection is the thickness of the slab and that is why it is not specified in the detector input, it is whatever the TissueInput tissue thickness is set to. This detector should be a function of Rho and Angle with angle being a polar angle as I described
above.
I specified Count=9 under Angle for TOfRhoAndAngleDetectorInput and set the slab thickness to be 10mm (set TissueInput, second Region, ZRange, Stop to be 10 and the third Region, ZRange, Start to be 10) within infile_one_layer_all_detectors.xml and
ran 50K photons. I also omitted the "axis image" from the matlab command for this plot. My results are shown here:
http://i.imgur.com/phyyW.jpg
A couple of notes:
1) When a photon goes out the slab beyond the last rho bin, we bin its contribution in the last bin. That is why there is a dark red band at the bottom left corner of the plot.
2) On the right hand half of the plot, there is no data. This is because there is air below the slab and so as the photon approaches the bottom surface of the slab from the inside, Fresnel will reflect the more grazing angles.
Hopefully you can reproduce these results.
Carole




Thanks so much for the reply. I am trying to duplicate the results but I am not successful. http://i.imgur.com/5x8YM.png
It seems like data from one Rho is repeated.
Can you please give me Matlab command u r using or any advice? I am setting my Delta=0.15707963267948966 and count=11, start=0, stop=1.5707963267948966.
I really do not know what is going on and I very much appreciate any help.




It sounds like you're specifying the TOfRhoAndAngleDetectorInput for angle correctly. So possibly you're having trouble with the Matlab commands to plot the results.
Similar to the steps dcuccia gave in his post above dated Wed at 12:14, the matlab commands are:
8) navigate Matlab to the "results" folder
9) Run the following:
h=fopen('TOfRhoAndAngle','rb');
transmittance=fread(h,[100,10],'double');
imagesc(log(transmittance));
Note that the sizes in the fread, [100,10] should match the number of rho and angle bins you've specified in TOfRhoAndAngleDetectorInput (one less than Count).




Thanks really so much. It worked out.
This is an amazing software with an outstanding support.
Alireza

