Save reflectance together with new asgn_float

229 views
Skip to first unread message

kris.dre...@gmail.com

unread,
Apr 24, 2024, 11:29:44 AM4/24/24
to mcx-users
Dear Qianqian,

Thank you for adding the long-awaited feature for the mua/mus/g/n 4xfloat format! :) 

We've implemented it on our side and executing mcx with "-b 1" and "-B aaaaaa" works for internal reflection combined with absorption at the volume borders.
However, if we try to include the detection of diffuse reflectance and photon exit positions, by adding "-w XV" and "-X 1", mcx runs as expected, but we get a shape mismatch when reading the mcx output. 

Did something change regarding the output shape of mcx when storing the reflectance or photon exit positions?

Best wishes and thank you!
Kris

Qianqian Fang

unread,
Apr 24, 2024, 12:19:36 PM4/24/24
to mcx-...@googlegroups.com, kris.dre...@gmail.com

hi Kris,

when you said "shape mismatch when reading the mcx output", were you referring to the shape of the detected photon data array when loading from the output file?

if so, what format did you ask mcx to save, and what function call did you use for loading the data?

if you use mcx in the command line, I suggest you use "-F jnii" because the output is a JSON file with explicit subfields for detected photons, making it easier to process.

if you use .mch, the loaded data must be properly parsed in order to understand which column of data correspond to what information you requested in -w flag.

for a simple test, I tried both mcxlab and mcx command line with -w XV and -X 1, the output look fine to me

./mcx --bench cube60 -w xv -X 1 -F jnii

cat cube60_detp.jdat | perl -pe 's/"[^"]{1000,}"/"..."/g'
{
    "MCXData":    {
        "Info":    {
            "Version":    1,
            "MediaNum":    2,
            "DetNum":    4,
            "ColumnNum":    6,
            "TotalPhoton":    1000000,
            "DetectedPhoton":    2984,
            "SavedPhoton":    2984,
            "LengthUnit":    1,
            "SeedByte":    0,
            "Normalizer":    200,
            "Repeat":    1,
            "SrcNum":    1,
            "SaveDetFlag":    48,
            "TotalSource":    1,
            "Media":    [{
                    "mua":    0,
                    "mus":    1.19209289550781e-07,
                    "g":    1,
                    "n":    1
                }, {
                    "mua":    0.00499999988824129,
                    "mus":    1,
                    "g":    0.00999999977648258,
                    "n":    1.3700000047683716
                }, {
                    "mua":    0.0020000000949949,
                    "mus":    5,
                    "g":    0.89999997615814209,
                    "n":    1
                }]
        },
        "PhotonData":    {
            "p":    {
                "_ArrayType_":    "single",
                "_ArraySize_":    [2984, 3],
                "_ArrayZipType_":    "zlib",
                "_ArrayZipSize_":    8952,
                "_ArrayZipData_":    "..."
            },
            "v":    {
                "_ArrayType_":    "single",
                "_ArraySize_":    [2984, 3],
                "_ArrayZipType_":    "zlib",
                "_ArrayZipSize_":    8952,
                "_ArrayZipData_":    "..."
            }
        }
    }
}


Qianqian

Kris --
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 on the web visit https://groups.google.com/d/msgid/mcx-users/9146c1fa-c805-4f46-9269-77b8269ade22n%40googlegroups.com.

kris.dre...@gmail.com

unread,
Apr 24, 2024, 12:39:35 PM4/24/24
to mcx-users
Hi Qianqian,

We tried it with both "-F jnii" and "-F bnii" which resulted in the output files with endings "_deip.jdat" and ".bnii". We are loading the files in Python with the jdata.load() function and for both options, it gives us a ValueError that looks like this: ValueError: cannot reshape array of size 85919032 into shape (64439274,)

Best wishes and Thank you!
Kris
 

Qianqian Fang

unread,
Apr 24, 2024, 4:19:27 PM4/24/24
to mcx-...@googlegroups.com, kris.dre...@gmail.com

can you update your jdata and bjdata python modules using the below command

python3 -m pip install --upgrade --force jdata bjdata

and rerun jd.load()? if it still give you an error, please upload your jdat or .bnii file in a shared drive such as dropbox or google drive and let me debug the issue.

kris.dre...@gmail.com

unread,
Apr 25, 2024, 5:05:48 AM4/25/24
to mcx-users
We are using the newest versions (bjdata==0.4.1 and jdata==0.6.0).

I uploaded the data to a Google Drive folder: https://drive.google.com/drive/folders/1K9LgFTss-zi2je4XoWRVpIjS7QF1nx9l?usp=sharing

Thank you for the help!

Qianqian Fang

unread,
Apr 25, 2024, 12:45:28 PM4/25/24
to mcx-...@googlegroups.com, kris.dre...@gmail.com

hi Kris,

I was able to reproduce the issue on my side.

it appears that you had added -Z 2 in the command line. this flag tells mcx to not apply any compression. because the output is JSON, a default base64 compression is required to convert binary to text, so this results in double base64 encoding, making the decoder confused.

I created a ticket for this bug

https://github.com/fangq/mcx/issues/219

and was able to fix it with the below patch

https://github.com/fangq/mcx/commit/ea67ea901c8be1e8bce5ce333bf133c9f84af31c

my recommendation is to use a -Z compression method that is not base64. if compression time is a concern, you can try lz4 (-Z 5) which is expected to be 20x faster than zlib.

the nightly build packages are being recompiled at this moment, please check out the nightly build folder later and test (or recompile mcx on your side). unfortunately github CI builds are failing right now because of ubuntu launchpad outatge: https://github.com/fangq/mcx/actions/runs/8836179945/job/24262240824

Qianqian

kris.dre...@gmail.com

unread,
Apr 26, 2024, 7:59:15 AM4/26/24
to mcx-users
Dear Qianqian, 
Thank you for fixing this, it worked for me now! :) 

The problem that I have now is that if I execute for example "mcx --bench cube60 -w xv -X 1 -F jnii", and load them with jdata.load(), I don't see any reflectance values. Usually, the reflectance values are negative, but there are no negative values in the output.

Best,
Kris

Qianqian Fang

unread,
Apr 26, 2024, 10:44:07 AM4/26/24
to mcx-...@googlegroups.com, kris.dre...@gmail.com

hi Kris,

saving dref (-X) requires the volume to have a layer of 0 padded right outside of the boundary where you want to receive dref readings. the built-in benchmark cube60 does not have this padded 0 layer, thus no dref output. you can manually add this 0 layer by using --json

mcx --bench cube60 -w xv -X 1 -F jnii --json '{"Shapes":[{"ZLayers":[1,2,0]}]}'

although you will also need to move up detector z positions by 1 voxel to avoid warnings.

Qianqian

kris.dre...@gmail.com

unread,
Jun 4, 2024, 9:36:03 AM6/4/24
to mcx-users
Hi Qianqian, we prepared a minimal working example that you can use to reproduce the behaviour I mentioned above.
There is also a notebook included that we use to visualize the results.


Thank you!

Best,
Kris

Qianqian Fang

unread,
Jun 4, 2024, 1:30:50 PM6/4/24
to mcx-...@googlegroups.com, kris.dre...@gmail.com

hi Kris,

as I mentioned in a recent thread posted by Seonyeong,

https://groups.google.com/g/mcx-users/c/GBDXxezzKEw/m/1v2gFsgBAwAJ

in non-label based media types, you can't set 0 to a voxel to indicate background, instead, you must set either mua or mus or both as NaN to indicate background medium.

because dref is only saved in the background voxels right next to non-zero voxels, you must pad a layer of NaNs in mua/mus float arrays in order for mcx to save dref data.

I loaded your vol.bin file, and I could not see 0 or NaN in the volume, so I assume that was the cause.

Qianqian

Liam Keegan

unread,
Jun 5, 2024, 10:26:48 AM6/5/24
to mcx-users
Hi Qianqian,

Thanks for looking into this issue, here is a notebook showing that the NaN padding layer is present in vol.bin and that with -b0 we get reflectance data, but when we use -b1 we don't get reflectance data:


Many thanks!

Liam

Qianqian Fang

unread,
Jun 5, 2024, 5:40:28 PM6/5/24
to mcx-...@googlegroups.com, Liam Keegan, kris.dre...@gmail.com

hi Liam and Kris,

I was able to load the volume again and can confirm that the NaN layer does exist.

To debug this issue, I recompiled mcx binary in the debug mode (make debug), and ran a single photon packet (-n 1) with either -b 0 or -b 1 with your command, here are the full logs of the first photon running in the two conditions (b1.log and b0.log).

I am not yet sure what was wrong, but I want to mention a few things while I was looking into this.

first thing is that when you set -b 0 to disable reflection, it not only disables the reflection at the boundary, but also throughout the entire domain, including between interior voxels. I see you have defined a structure for n per voxel, optically, this may be a very different problem when you model light reflections vs ignoring it between spatial heterogeneous voxels. but I am not certain this is the case here. I set the refractive index to 1.3 constant across the domain, I also saw very different reflection behavior.

if you do a diff to the two log files, the first 186 lines are the same, it starts to differ when the photon moves from voxel [76 29 19] to [76 29 20], where the per-voxel defined refractive index, n, is different in the two voxels. Therefore, in the -b 1 case, reflection/refraction calculations are invoked while the other one does not. From what I can tell, the -b 1 case, most photons exits from other surfaces but in -b 0, most of those are back reflected.

again, I do not fully understand what has caused this difference - whether this difference is true when considering or ignoring n-mismatch, or it is a bug in the code. I just want to share this perhaps you could also read the log, or use the debug mode to investigate this in parallel. you could also plot the trajectories and see why photons are mostly forward-scattered in -b 1 but not in -b 0. I suspect that there might be something related to the boundaries of the low mua/low mus voxels (mua=1e-5, mus=1e-5).

Qianqian

b1.log
b0.log

Liam Keegan

unread,
Jun 27, 2024, 9:37:23 AM6/27/24
to mcx-users
Hi Qianqian,

Sorry for the delay in responding - I now have a simple test case similar to your example, with a uniform volume with n=1.3 for all voxels and for the background, where a single photon exits at the z=0 surface padded with NaNs.

With b=0 the photon exits voxel [93 12 0] - and I can read the resulting reflectance.

With b=1 the photon is instead reflected, with the message `ref total ref=-inf`.

Looking into this it turns out that both g and n are zero at these padded voxels - I think this is because g/n are not read if mua/mus are nan: https://github.com/fangq/mcx/blob/f959c71f6d0ba6471a15a6fd016480fd4b5efd29/src/mcx_utils.c#L3526-L3538

I tried to modify this so that g/n are read: https://github.com/fangq/mcx/pull/224

With this change the single photon example now has the same behaviour with b=1 and b=0, however when I tried a full simulation I produced many errors of the form:

ERROR: should never happen! mediaid=0 idx1d=1275 isreflect=1 gcfg->doreflect=1 n1=1.299805 n2=1.299805 isdet=0 flipdir[3]=2 p=(45.095932 39.130569 1.000000)[45 39 0]

So I'm not sure if my change is helpful or not.

Many thanks for your help!

Best wishes,

Liam

Liam Keegan

unread,
Jun 27, 2024, 9:45:44 AM6/27/24
to mcx-users
I couldn't upload the logs here so I created an issue with the log files: https://github.com/fangq/mcx/issues/225

Sorry for the multiple messages.

Best wishes,

Liam

Reply all
Reply to author
Forward
0 new messages