For MEDIA_ASGN_F2H where mua/mus are nan ensure that the provided g/n…

[lkeegan]

… values are nontheless read

Check List

Before you submit your pull-request, please verify and check all below items

• You have only modified the minimum number of lines that are necessary for the update; you should never add/remove white spaces to other lines that are irrelevant to the desired change.
• You have run make pretty (requires astyle in the command line) under the src/ folder and formatted your C/C++/CUDA source codes before every commit; similarly, you should run python3 -m black *.py (pip install black first) to reformat all modified Python codes, or run mh_style --fix . (pip install miss-hit first) at the top-folder to format all MATLAB scripts.
• Add sufficient in-code comments following the doxygen C format
• In addition to source code changes, you should also update the documentation (README.md, mcx_utils.c and/or mcxlab.m) if any command line flag was added or changed.

If your commits included in this PR contain changes that did not follow the above guidelines, you are strongly recommended to create a clean patch using git rebase and git cherry-pick to prevent in-compliant history from appearing in the upstream code.

Moreover, you are highly recommended to

• Add a test in the mcx/test/testmcx.sh script, following existing examples, to test the newly added feature; or add a MATLAB script under mcxlab/examples to gives examples of the desired outputs
• MCX's simulation speed is currently limited by the number of GPU registers. In your change, please consider minimizing the use of registers by reusing existing ones. CUDA compilers may not be able to optimize register counts, thus require manual optimization.
• Please create a github Issue first with detailed descriptions of the problem, testing script and proposed fixes, and link the issue with this pull-request

Please copy/paste the corresponding Issue's URL after the below dash

[fangq]

@lkeegan, can you test my above fix in 80b5794 ?

this should allow mcx to read g/n regardless mua/mus=nan or not.

[lkeegan]

With 80b5794 I see very different behaviour in the MEDIA_ASGN_F2H case with boundary padding.

I think this is because

• the nan valued mua/mus are set to zero:

mcx/src/mcx_utils.c

Lines 3526 to 3532 in 80b5794

	if (f2h[0] != f2h[0] || f2h[1] != f2h[1]) { /*if one of mua/mus is nan in continuous medium, convert to 0-voxel*/
	cfg->vol[i] = 0;
	if (cfg->mediabyte == MEDIA_AS_F2H) {
	continue;
	}
	}

• but then they are overwritten with their original (nan) values:

mcx/src/mcx_utils.c

Lines 3534 to 3539 in 80b5794

	if (cfg->mediabyte == MEDIA_ASGN_F2H) {
	cfg->vol[i] = mcx_float2half2(f2h);
	f2h[0] = val[(i << 2) + 2]; // g
	f2h[1] = val[(i << 2) + 3]; // n
	cfg->vol[i + datalen] = mcx_float2half2(f2h);
	} else {

So I think the g/n values are now being read, but the boundary mua/mus values are now nan instead of 0

[fangq]

sorry for taking a while to get back to this - let me know if the above commit fixes the issue.

[lkeegan]

Thanks! With this commit the behaviour is the same as #224 - the ref total ref=-inf issue is fixed, but I still don't get any reflectance with b=1:

https://gist.github.com/lkeegan/ee77309c2b48883278fd9716c95d1c34

And there are also many error messages of the form:

ERROR: should never happen! mediaid=0 idx1d=88F isreflect=1 gcfg->doreflect=1 n1=1.299805 n2=1.299805 isdet=0 flipdir[3]=2 p=(31.562723 18.700874 1.000000)[31 18 0]

(see #225, https://groups.google.com/g/mcx-users/c/bToGluYYdao/m/5XGk9jeHAAAJ)

[fangq]

ok, let me do more test on this.

[kdreher]

I was wondering if there are any updates on this.

[fangq]

@lkeegan and @kdreher, sorry for my late response.

I could not find the files vol_bg.bin and vol_bg.json in your jupyter notebook, so I wrote a simple test case based on a built-in demo

clear cfg;
cfg.nphoton = 1e7;
cfg.vol = uint8(ones(60, 60, 60));
cfg.srcpos = [30 30 1];
cfg.srcdir = [0 0 1];
cfg.gpuid = 1;
% cfg.gpuid='11'; % use two GPUs together
cfg.autopilot = 1;
cfg.prop = [0 0 1 1; 0.005 1 0 1.37];
cfg.tstart = 0;
cfg.tend = 5e-9;
cfg.tstep = 5e-9;

mua1 = 0.003;
mua2 = 0.1;
mua = single(reshape(repmat([mua1:(mua2 - mua1) / (60 - 1):mua2]', 60, 60), 60, 60, 60));

mus1 = 1.0;
mus2 = 5.0;
mus = single(reshape(repmat([mus1:(mus2 - mus1) / (60 - 1):mus2]', 60, 60), 60, 60, 60));

cfg.vol = reshape([mua(:)'; mus(:)'], [2 60 60 60]);
cfg.vol(3, :, :, :) = 0;
cfg.vol(4, :, :, :) = 1.37;

cfg.vol(1,:,:,1)=nan;
cfg.issaveref=1;
cfg.isreflect=1;

flux = mcxlab(cfg);
mcxplotvol(log10(flux.dref))

with the above script, I was able to get non-empty dref, see screenshot below - can you reproduce this in your environments (mcx binary+python)?

[fangq]

closing for now as the issue can no longer be reproduced.

feel free to reopen if the problem persists.

[lkeegan]

Thanks, I'll try to reproduce your example script with python+mcx.

To reproduce the original reported issue, here is a zipfile with vol_bg.bin and vol_bg.json:

vol_bin_and_json.zip

And the command line args for mcx:

mcx -f vol_bg.json -O F -a 1 -F jnii -Z 2 -b 1 --saveref 1 -s b0_bg

Running the above results in no reflectance data and many error messages of the form:

ERROR: should never happen! mediaid=0 idx1d=1981 isreflect=1 gcfg->doreflect=1 n1=1.299805 n2=1.299805 isdet=0 flipdir[3]=2 p=(49.131416 54.812988 1.000000)[49 54 0]

If I set -b 0 instead then it runs without error messages and produces the expected reflectance data.

Many thanks for your help!