Converting between 'energy' and 'flux' outputs

[Jessica Lund]

Hi Dr. Fang,

I'm working on a look-up table and I need both the energy distribution and the flux distribution in my volume. 

I've seen in a few other posts here that Energy = Flux*ua*unitinmm so I've been doing one simulation for each set of optical properties and using that conversion factor to get both outputs.

However, I noticed that when ua<0.004 mm^-1, the 'flux' output drops. Should I assume that the 'energy' output is correct and use that to get the flux distribution?

Thank you!

[Jessica Lund]

Also, I noticed that it is dependent on the voxel size. In my original post, unitinmm = 0.25, and this is with unitinmm = 0.2 

So it seems that the cutoff = (1/unitinmm)/1000

[Fang, Qianqian]

Hi Jessica,

Thanks for reporting this. Can share your script how you created this plot?

The cutoff mua=0.001 is hardcoded here

https://github.com/fangq/mcx/blob/v2024.2/src/mcx_core.cu#L2104C42-L2104C48

As the result of a commit based on a patch proposed by Ernesto Pini

https://github.com/fangq/mcx/pull/164

After the simplified patch (https://github.com/fangq/mcx/commit/0b9884334efcfdce13c24e83b2ece988166535e0), the test script in the above pull request produces smooth transition from high to low mua values (which I just confirmed again).

I am now curious why your plot is not continuous - can you provide your script how you created this plot?

Qianqian

---

From: mcx-...@googlegroups.com <mcx-...@googlegroups.com> on behalf of Jessica Lund <jessic...@gmail.com>
Sent: Friday, November 15, 2024 12:13 PM
To: mcx-users <mcx-...@googlegroups.com>
Subject: [mcx-users] Re: Converting between 'energy' and 'flux' outputs

 

--
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 view this discussion visit https://groups.google.com/d/msgid/mcx-users/b1dc4404-32b4-40b8-ba6f-407850e79d97n%40googlegroups.com.

[Jessica Lund]

Here's my script:

clear; clc;

n = 1.4;
g = 0.9;

% length of voxel in mm
cfg.unitinmm = 0.25;

nR = 20/cfg.unitinmm;
nz = 7/cfg.unitinmm;

% volume
cfg.vol = uint8(ones(nR,nR,nz));

% source parameters
probe_NA = 0.22;
probe_r = 0.5; % mm
cfg.srctype = 'disk';
cfg.srcpos = [nR/2 nR/2 nz];
cfg.srcparam1 = [probe_r/cfg.unitinmm 0 0 0];
cfg.srcparam2 = [0 0 0 0];
f = -probe_r/tan(asin(probe_NA/n))./cfg.unitinmm; % focal length
cfg.srcdir = [0 0 -1 f];
cfg.issrcfrom0 = 1;

% other simulation inputs
cfg.prop = [5 0 1 1.5; 0 0 g n];
cfg.nphoton = 1e7;
cfg.bc = 'rrrrra';
cfg.savedetflag = 'dpxvw';
cfg.gpuid = 1;
cfg.autopilot = 1;
cfg.tstart = 0;
cfg.tend = 5e-9;
cfg.tstep = 5e-9;
cfg.debuglevel = 'P';
cfg.isnormalized = 0;

% start the comparison
num_points = 30;

% choose opt props
min_ua = 0.0025;
max_ua = 0.0055;
ua = linspace(min_ua,max_ua,num_points);
us = 25;

Freal = zeros(1,num_points);
Fcheck = zeros(1,num_points);
Ereal = zeros(1,num_points);
Echeck = zeros(1,num_points);

h = waitbar(0, 'Starting');
cfg.prop(2,2) = us;
for i = 1:num_points
    cfg.outputtype = 'flux';
    cfg.prop(2,1) = ua(i);
    [flux_em, ~, ~, ~] = mcxlab(cfg);
    Freal(i) = sum(flux_em.data./cfg.nphoton,'all');
    Echeck(i) = sum(flux_em.data./cfg.nphoton.*ua(i).*cfg.unitinmm,'all');
   
    cfg.outputtype = 'energy';
    [energy_em, ~, ~, ~] = mcxlab(cfg);
    Ereal(i) = sum(energy_em.data./cfg.nphoton,'all');
    Fcheck(i) = sum(energy_em.data./cfg.nphoton./(ua(i).*cfg.unitinmm),'all');
    waitbar(i/num_points,h, sprintf('Percent completed: %0.3f %%', i/num_points*100));
end    
delete(h)
%%
figure; plot(ua,Ereal); ylabel('energy'); xlabel('ua (mm^{-1})');
hold on; plot(ua,Echeck,LineStyle='--');
legend("outputtype='energy'", "outputtype='flux'")
grid on
figure; plot(ua,Freal); ylabel('flux'); xlabel('ua (mm^{-1})');
hold on; plot(ua,Fcheck,LineStyle='--');
legend("outputtype='flux'", "outputtype='energy'")
grid on

[Fang, Qianqian]

Hi Jessica,

Thanks for the test script. I see the problem.

I believe this is a result of a patch previously introduced in https://github.com/fangq/mcx/pull/164, specifically,

https://github.com/fangq/mcx/commit/481c4b5e5ca10e3b0b46b5a42febb8a7399355fe

It appears that the threshold (0.001/unitinmm) to apply the low-mua asymptotic formula is too high, creating an error for not-so-low mua values.

I eventually lower this threshold to EPS (~1e-7) by changing to the line to

 weight = (prop.mua < EPS) ? (w0 * len) : __fdividef(w0 - p.w, prop.mua);

If I run your script and lower the mua range to bracket ~1e-7, you can only see the discrepancy bellow 1e-6.

Ernesto: I reran your test script for #164, now even with the prop.mua<EPS setting, I was able to get a smooth asymptotic line - I am now wondering what was the thought for setting the threshold to 0.001?

Qianqian

---

From: mcx-...@googlegroups.com <mcx-...@googlegroups.com> on behalf of Jessica Lund <jessic...@gmail.com>

Sent: Saturday, November 16, 2024 11:33 AM
To: mcx-users <mcx-...@googlegroups.com>
Subject: Re: [mcx-users] Re: Converting between 'energy' and 'flux' outputs

 

To view this discussion visit https://groups.google.com/d/msgid/mcx-users/0928f1ff-af9b-43a3-966b-78890ad0f6ddn%40googlegroups.com.

[Fang, Qianqian]

Jessica, I also want to point out that the following mcxlab call you used

[flux_em, ~, ~, ~] = mcxlab(cfg);

Can lead to inefficient computation. Please see my reply point#2 in this post: https://groups.google.com/g/mcx-users/c/Bk75GmKug-Q/m/VjjGtpwwBQAJ

Although matlab supports ~, inside called function, there is no way to know which output should be ignored. mcxlab simply uses nargout to decide what output to bypass. But using ~ won't shorten nargout (in the above call, nargout is 4).

If you call mcxlab using your above statement, mcxlab will have to compute all 4 outputs, but you only need the first output. This can slow down your computation. Use the below statement instead

flux_em = mcxlab(cfg);

---

From: 'Fang, Qianqian' via mcx-users <mcx-...@googlegroups.com>
Sent: Monday, November 18, 2024 11:46 PM
To: mcx-users <mcx-...@googlegroups.com>
Cc: Ernesto Pini <pi...@lens.unifi.it>

To view this discussion visit https://groups.google.com/d/msgid/mcx-users/BYAPR06MB41032715F858E4114DF569B7E4202%40BYAPR06MB4103.namprd06.prod.outlook.com.

[Ernesto Pini]

Dear Prof Fang,

Sorry for the late reply,

The idea of the cutoff was to set it when I observed a deviation in the fluence from the nominal value.
Originally the cutoff was set at mua/mus < 0.001 since scattering also plays a role in this, but I see that you changed to mua < EPS to account for mus=0.
Is that correct?
The only problem I can think of is that now the cutoff is absolute and not relative (adimensional).

Kind regards
Ernesto

[Jessica Lund]

Thank you!