Normalization of Flux in MCX

[Timothy Sowers]

Hi Dr. Fang, 

I'm running simulations in MCXLAB in which I'm trying to compare flux at different depths for a photoacoustic application.  Based on the descriptions in the MCX Readme, the flux is already normalized such that the source has a value of 1.  

I wanted to double check my understanding of this with some basic calculations using the flux distributions I'm getting from my simulations, but they are not working out.  

My simulations space is a 30 x 30 x 30 mm volume, and I'm using a voxel size of 0.1 mm.  My source is a planar source that is 1.2 x 16 mm in size.  I'm checking the flux using the following approach.  Calculate the average flux across all the elements directly in front (the first element only) of the source, then multiply by the area (1.2*16) of the planar source.  Based on the description in the readme, I'd think this value should be 1.  However, the value I'm getting is 1.87.  A number slightly lower than 1 would make sense to me due to losses from reflection or backwards scattering (and because the direction of my source is sideways by ~50 degrees), but this value seems wrong to me

Is there a physical phenomena I'm not considering or is my logic somehow not correct here?  

I've attached the MATLAB code I've used to run the simulation (simulation.m), an image of a 2D slice of the volume directly in front of the source (image.png), the short loop I used to calculate the value of 1.87 (calc.m), and the output of MCXLAB before running is also pasted below.  

Thanks, 

Tim

Launching MCXLAB - Monte Carlo eXtreme for MATLAB & GNU Octave ...

Running simulations for configuration #1 ...

mcx.nphoton=1e+08;

mcx.dim=[300 300 300];

mcx.unitinmm=0.1;

mcx.srctype='planar';

mcx.srcpos=[239 70 1];

mcx.srcparam1=[12 0 0 0];

mcx.srcparam2=[0 160 0 0];

mcx.srcdir=[-0.438371 0 0.898794];

mcx.gpuid=1;

mcx.autopilot=1;

mcx.medianum=3;

mcx.tstart=0;

mcx.tend=5e-09;

mcx.tstep=5e-09;

###############################################################################

#                      Monte Carlo eXtreme (MCX) -- CUDA                      #

#          Copyright (c) 2009-2016 Qianqian Fang <q.fang at neu.edu>          #

#                             http://mcx.space/                               #

#                                                                             #

#         Computational Imaging Laboratory (CIL) [http://fanglab.org]         #

#            Department of Bioengineering, Northeastern University            #

###############################################################################

#    The MCX Project is funded by the NIH/NIGMS under grant R01-GM114365      #

###############################################################################

$Rev::       $ Last $Date::                       $ by $Author::              $

###############################################################################

- variant name: [Fermi] compiled for GPU Capability [100] with CUDA [7000]

- compiled with: RNG [POSIX erand48] with Seed Length [4]

- this version CAN save photons at the detectors

GPU=1 (GeForce GT 730) threadph=24414 extra=256 np=100000000 nthread=4096 maxgate=1 repetition=1

initializing streams ... init complete : 14 ms

requesting 3072 bytes of shared memory

lauching MCX simulation for time window [0.00e+00ns 5.00e+00ns] ...

simulation run# 1 ... kernel complete:   235123 ms

retrieving fields ... transfer complete: 235193 ms

data normalization complete : 235194 ms

normalizing raw data ... normalization factor alpha=199.999985

simulated 100000000 photons (100000000) with 4096 threads (repeat x1)

MCX simulation speed: 425.34 photon/ms

total simulated energy: 100000000.00 absorbed: 25.93493%

(loss due to initial specular reflection is excluded in the total)

flux1100 = 

    data: [300x300x300 single]

    stat: [1x1 struct]

[Timothy Sowers]

Hi Dr. Fang, 

As a quick follow up to this, I have another image that makes the point in a way that is more straightforward.  It's attached (image2).  You can see that the max flux right in front of the source is about 0.12.  If the source had a normalized value of 1, then the flux should be 1/(1.2 mm x 16 mm) = 0.052 mm^-2, correct?  Note that I've already multiplied the flux.data by cfg.tend in this image to convert to fluence (I only used one step).  

Any insight would be appreciated.  

Thanks, 

Tim

[Timothy Sowers]

I forgot to attach

[Timothy Sowers]

Hi Dr. Fang, 

If you get the time I'd appreciate any insight you can give into this problem.  I've been comparing different designs for a photoacoustic device.  Because of this, I've been able to compare relative fluence distributions between different designs, and not worry about the absolute values.  However, for future projects the absolute fluence values are great to have.  

So far, the simulations comparing different designs are proving very informative and non-intuitive.  Thanks for putting the program together.  It's quite powerful and useful.  

Tim

[Qianqian Fang]

On 03/16/2017 10:37 AM, Timothy Sowers wrote:

Hi Dr. Fang, 

If you get the time I'd appreciate any insight you can give into this problem.  I've been comparing different designs for a photoacoustic device.  Because of this, I've been able to compare relative fluence distributions between different designs, and not worry about the absolute values.  However, for future projects the absolute fluence values are great to have.  

So far, the simulations comparing different designs are proving very informative and non-intuitive.  Thanks for putting the program together.  It's quite powerful and useful. 

hi Tim

I apologize for the delay in getting back to you.

I first would like to make some clarifications on the quantities resulted from
mcx simulations.

The term "flux" has many meanings in physics and is indeed very confusing
(I am sorry for choosing that ambiguous name in the documentation).

What MCX/MMC really outputs is the "fluence rate Phi(r,t)" (i.e. the first output of mcxlab.m).
It has a unit of "W/m^2", referring to the expected number of photons, or
Joule of energy passing by position r per unit time at time t, regardless of direction.
This is also the TPSF generated at each voxel.

If you multiply the fluence rate by the time-gate width (cfg.tstep) and sum all
the time gates, that yields "fluence" - Phi(r) (i.e. the so called CW fluence),
which has a unit of J/m^2.

MCX can also output "current density" (J(r,t), unit W/m^2, same as Phi(r,t)) -
referring to the expected number of photons or Joule of energy flowing
through a unit area pointing towards a particular direction per unit time.
The current density can be calculated at the boundary of the domain by two means:

1. using the detected photon partial path output (i.e. the second output of mcxlab.m),
one can compute the total energy E received by a detector, then one can
divide E by the area/aperture of the detector to obtain the J(r) at a detector
(E should be calculated as a function of t by using the time-of-fly of detected
photons, the E(t)/A gives J(r,t); if you integrate all time gates, the total E/A
gives the current I(r), instead of the current density).

2. use the new -X 1 or --saveref/cfg.issaveref option in mcx to enable the
diffuse reflectance recordings on the boundary. the diffuse reflectance
is represented by the current density J(r) flowing outward from the domain.

The current density has, as mentioned, the same unit as fluence rate,
but the difference is that J(r,t) is a vector, and Phi(r,t) is a scalar. Both measuring
the energy flow across a small area (the are has direction in the case of J) per unit
time.

You can find more rigorous definitions of these quantities in Lihong Wang's
Biomedical Optics book, Chapter 5.

Here comes the confusion, because the current density sometimes is
also called radiative flux or flux (unit W/m^2) in some fields of physics,
so I labeled the output of mcx as flux in the documentation, but it was
indeed an incorrect choice of terminology - a flux is typically a vector,
and a fluence/fluence rate is a scalar.

In short, the first output of mcxlab is the fluence rate Phi(r,t), the second
output detps can derive the current density J(r,t) at the detectors.

After normalization, the fluence rate output should match the fluence
rate derived from the RTE (in the case of large photons), and should also
match nicely with the diffusion approximation in the diffusive regime
when the DA is valid (i.e. high scattering). So, my effort of validation
focuses primarily on comparing the MCX fluence/fluence rate output
with the diffusion equation's Green's function (assuming unitary energy
source) in the diffusive regime.

having these output definitions in mind, I think the observations you made
seem to make sense to me.

if your medium has no scattering (like air or water), then, if you launch
1 Joule of photons at your source, any area/voxel in front of the source
should receive exactly 1 J / area fluence, because a photon only passes by
a voxel one-way one time, and never returns.

However, for a scattering medium, a photon can pass by a voxel wall
many times due to scattering. Thus, the total net energy flowing through
any voxel boundary is more than what it is supposed to receive when
there is no scattering.

I will update the documentation to reflect the more rigorous definitions of
these quantities. In the meantime, if this clarification does not help
understanding the observations you made, let's continue this discussion.

Qianqian

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To post to this group, send email to mcx-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mcx-users.
For more options, visit https://groups.google.com/d/optout.

[Jiamengyu]

hi Dr. qianqian, 

thanks for your specific explanation about the output of MCX! The term of 'flux' in the document also made me confused before.

Sincerely yours,

Mengyu Jia

Tianjin university, China

发自我的 iPhone

[Timothy Sowers]

Hi Dr. Fang, 

Thanks for the quick response.  It was complete and very clear.  

Your explanation makes perfect sense - that in a scattering medium the current density (which I'm calculating) would not be equivalent to the fluence rate (which mcxlab gives in flux.data).  

Over the next few days, I'll run some simulations on non-scattering media and also put detectors into my current simulations to double check my understanding.  I don't think further discussion will be needed, although if I run into anything unanticipated in the simulations or come up with any other questions I'll be sure to ask.  

Thanks for the help!  

Tim