Strange artifact in MCXLab

[Lukai Wang]

Dear Prof. Fang:

Hope this letter find you well! I am a beginner in Monte Carlo simulation. When I was using MCXLab (v2023) in MATLAB, I noticed a weird artifact. I was simulating a pencil beam perpendicularly incident into a homogeneous medium (0.001 0.04 0.9 1.37) in a 2D simulation modified from your example code demo_matlab_2d.m with the following code and I encountered such artifacts:

%==========================================================================

% A sample 2D MCX simulation

%

% Author: Qianqian Fang <q.fang at neu.edu>

%==========================================================================

% you must define a 3D array with one singleton dimension (with length 1).

% unfortunately if you define z as singleton, matlab will make the array 2D

% instead of 3D, so we have to purmute it to make the 1st dimension

% singleton.

% only clear cfg to avoid accidentally clearing other useful data

clear cfg

cfg.nphoton=1e6;

% Test the situation with only on medium

cfg.vol=permute(uint8(ones(100,100)), [3,1,2]); % from 2d to 3d

% cfg.vol(1,30:70, 10:50)=2;

cfg.issrcfrom0=1;

cfg.srctype='pencil';

cfg.srcpos=[0,50,0]; % src position must be located in the 2D plane

cfg.srcdir=[0 0 1]; % src dir must align in the plan (y-z in this case)

cfg.gpuid=1;

cfg.autopilot=1;

cfg.gscatter=100;

myprop=[0.001 0.04 0.9 1.37];

cfg.prop=[0 0 1 1; myprop];

cfg.tstart=0;

cfg.tend=1e-8;

cfg.tstep=1e-8;

flux=mcxlabcl(cfg);

cw=squeeze(sum(flux.data,4));

figure;

imagesc((abs(cw)))

figure;

plot(abs(cw(51,:)));

The above two figures are the 2D fluence distribution and 1D fluence distribution

along the pencil beam direction. As we can see, there is an obvious oscillation of

fluence with respect to the coordinate.

In this simulation, the voxel size is around 4% of the mean free path and 0.4%

transport mean free path. And I found out that increasing voxel size will denoise

this artifact. So is this artifact caused by the small voxel size? Or is there other reasons?

Thank you for your time. Your kindly help would be deeply appreciated.

Best,

Lukai

[Qianqian Fang]

hi Lukai,

it seems you were not using mcxlab, but the OpenCL version mcxlabcl.

unfortunately I am not able to reproduce this issue. I tried both Linux and Windows, and the output looks fine to me.

see attached screenshot: 

- left: windows, 

- right: Linux, 

    - top-row on the right: running on the CPU (AMD Ryzen7 4800H with pocl runtime), 

    - bottom-row: on the NVIDIA GPU (RTX2060).

please report your OS and OpenCL device - use mcxlabcl('gpuinfo') to print and use cfg.gpuid to choose.

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 view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/53cf99ca-6721-4832-af7c-952ea44ab095n%40googlegroups.com.

[Lukai Wang]

Dear Prof. Fang:

Thank you for the detailed response. I was using V2020 yesterday. Today I have made the following attempts with two versions of MCX (2023 and 2020) and both mcxlab (calls mcxcl.m) and mcxlabcl on my office computer (Windows) and laptop (MacOS). mcxlab and mcxlabcl always give the same result. Seems like it's openCL device that matters. Only Intel(R) Core(TM) i5-1038NG7 CPU in my Macbook seems to work. In addition, the latest version V2023 seems to give an incorrect sudden rise in flux. 

1. Windows - Intel(R) HD Graphics 530:

V.2023 / mcxlab: No artifact, but a very uncommon raise appears at around 30, which is obviously not correct. mcxlabcl gives same results.

V.2020 / mcxlab: Significant artifact. mcxlabcl gives same results.

2. MacOS - Intel(R) Iris(TM) Plus Graphics / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz

V.2023 / mcxlab / Intel(R) Iris(TM) Plus Graphics:  Uncommon raise in flux.

V.2023 / mcxlab / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz: Seems correct

V.2020 / mcxlab / Intel(R) Iris(TM) Plus Graphics: Artifact 

V.2020 / mcxlab / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz: Seems correct.

So I believe the reason for the artifact is the OpenCL device. 

Although the result by V.2020 / mcxlab / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz is continuous and smooth. I still have some doubt, I changed anisotropy g from 0 to 0.9 in the code and the fluences seem almost the same:

But when I was comparing with Prof. Lihong Wang's code MCML (https://omlc.org/software/mc/mcml/index.html) with the same parameters, the maximum of flux appears at the depth of transport mean free path instead of the boundary and g plays an important role as follows:

Could you give me some advice on this difference as well? Your kindly help would be greatly appreciated.

Best,

Lukai

[Qianqian Fang]

hi Lukai,

two comments:

1. stop using Intel HD/Iris GPUs with mcxcl - I have noticed that Intel HD/Iris GPUs had failed a subset of my built-in tests (https://github.com/fangq/mcxcl/blob/master/test/testmcx.sh), I believe these issues were caused by Intel GPU's OpenCL runtime, but haven't been able to pin-point what was the problem - all other CPUs/GPUs we have tested mcxcl have all passed, except HD/Iris integrated GPUs.

my group also purchased a discrete Intel GPU A770 earlier this year, and it also runs fine for all tests of mcxcl. so the issue seems to be specific to HD/Iris iGPUs. at this moment, please just use mcxcl on your CPU (feel free to run the testmcx.sh script in a bash shell to verify the GPU is working).

2. please do not use 2D simulations for any quantitative studies - the 2D simulation is not physically inspired - it is a faster simulation mode that I used for testing purposes (and training data generation for some of my machine-learning based works), but it is not how photons propagates in space. although MCML outputs 2D fluence output, it still simulates photon's movement in full 3D space. you should use a 3D volume with mcxcl if you want to match mcml output.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/179baf99-4838-45fc-beb3-1622a9a100c1n%40googlegroups.com.

[Qianqian Fang]

another quick comment - I know many users attempted to run mcxcl/mmc with the integrated GPU because it is abundantly available - what I want to say is it is not the best hardware to take advantage of the scalability of mcxcl/mmc. even it can run correctly, its speed up over the CPU is still quite small, and mcxcl/mmc can be a lot faster (and more reliable) if you give it a discrete GPU.

if you want to help debugging this, can you set cfg.optlevel=0 in your script and see if the artifact is still there on your iGPUs? mcxcl has 4 optimization levels (see our 2018 paper), and currently the default optimization level is 1 (which enables -cl-mad-enable flag and use the native_* versions of the math functions) - I don't know if these were not implemented properly in Intel's HD/Iris runtimes, but disabling all optimization by setting cfg.optlevel=0 could tell.

Qianqian

[Lukai Wang]

Thank you Prof. Fang. I modified the previous 2D test code into the following 3D code and tested with g = 0 / 0.9 , but the results are still almost the same:

% only clear cfg to avoid accidentally clearing other useful data
clear cfg

cfg.nphoton=1e6;

% Test the situation with only one medium
cfg.vol=uint8(ones(20,20,100)); % from 2d to 3d
cfg.issrcfrom0=1;
cfg.srctype='pencil';

cfg.srcpos=[10,10,0];   % src position

cfg.srcdir=[0 0 1];    % src dir must

cfg.gpuid=2;
cfg.autopilot=1;
cfg.gscatter=100;

myprop=[0.001 0.04 0 1.37];

cfg.prop=[0 0 1 1; myprop];
cfg.tstart=0;
cfg.tend=1e-8;
cfg.tstep=1e-8;

flux = mcxlab(cfg, 'opencl');
% flux=mcxlabcl(cfg);


cw=squeeze(sum(flux.data,4));
figure;
imagesc((abs(squeeze(cw(11,:,:)))));
title('flux, g = 0.9', 'fontsize', 20);

figure;
plot(abs(squeeze(cw(11,11,:))));
title('flux along the center, g = 0.9', 'fontsize', 20);

Is there something wrong with my setup? Your kindly help is really appreciated.

In addition, I tried cfg.optlevel=0 in my script with my iGPUs, but the artifacts remain the same. 

Best regards,

Lukai

[Qianqian Fang]

are you sure your MCX and MCML simulations used the same setting?

for example, what is the domain size/length scale of your MCML set up? does mcx volume match that?

for low-scattering cases like this one, the plot should really just follow the Beer-Lambert law and it should be fairly easy to verify

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/6288b745-1650-40dd-a71f-899bbe416301n%40googlegroups.com.

[Lukai Wang]

Yes I used the same setting:

The domain size is 10mm x 10mm x 100mm, with voxel size of 1mm. 

The medium parameter is: μs = 0.04mm^-1, μa = 1e-3mm^-1, g = 0, n = 1.37. 

In MCX I get an exponential distribution which fits the Beer Lambert law, while I get a distribution with a maximum at mean free path in MCML:

MCX:

MCML:

I am confused why the results are different with the same setting. I think a peak at around transport mean free path makes sense because of backscattering in tissue.

Best,

Lukai

[Qianqian Fang]

what is the y-axis of your mcml plot? the two curves have totally different shapes, making me think that you might have plotted different quantities.

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/edc458bb-815e-462b-b6c8-1be031200508n%40googlegroups.com.

[Lukai Wang]

It is normalized fluence.

[Qianqian Fang]

I still don't understand the rationale why you expect the fluence to have a peak at a depth of mean-free path, I also don't have an explanation why your mcml output did not show an exponential decay ... you will have to understand what you were simulating with mcml and build an equivalent one with mcx. right now, I don't think they are simulating the same thing

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/fc5038ca-0613-400c-bcae-abe4f65bc6bcn%40googlegroups.com.

[Lukai Wang]

Thank you very much Prof. Fang. I will check more carefully. Right now I think I am simulating the fluence of a pencil beam propagating in a homogeneous medium with mcml and mcx. 

My understanding of a peak of fluence at a depth of mean-free path is that photons get scattered through the propagation and the peaks appears around the depth of mean-free path, which is shown in Chapter 3 of 《Biomedical Optics: Principles and Imaging》 by Prof. Lihong Wang (http://www.morrisvanegas.com/ref/eece5648/%5BLihong_V._Wang,_HsinI_Wu(auth.)%5D_Biomedical_Opti(BookZZ.org).pdf) as follows: 

[Qianqian Fang]

hi Lukai,

I suspected that's where it came from.

I saw a couple of problems:

1. the peak at 1/mus'=l'tr is a result of diffusion due to high scattering. however, your benchmark is non-diffusive where mus is extremely low. because of that, photons are largely ballistic and the dominant effect is the exponential decay along the trajectories.

2. the depth-resolved fluence in your cited reference is not just the fluence along the incident line-path of the pencil beam, rather, it is the integration along the radial direction (in a voxelated domain, it is the summation along x/y planes) - see Eq. 4.20 in MCML's manual. It is mean to be an approximation to the 1D decay profile when you have an infinite plane-wave passing through a semi-infinite slab.

please consider modifying your code to see this effect by summing the x/y dimensions:

% only clear cfg to avoid accidentally clearing other useful data
clear cfg

cfg.nphoton=1e6;
% Test the situation with only one medium


cfg.vol=uint8(ones(40,40,100)); % from 2d to 3d
cfg.issrcfrom0=1;
cfg.srctype='pencil';

cfg.srcpos=[20,20,0];   % src position


cfg.srcdir=[0 0 1];    % src dir must



cfg.gpuid=1;
cfg.autopilot=1;
cfg.unitinmm=0.1;

myprop=[0.01 10 0.9 1.37];



cfg.prop=[0 0 1 1; myprop];
cfg.tstart=0;
cfg.tend=1e-8;
cfg.tstep=1e-8;


flux = mcxlab(cfg);
% flux=mcxlabcl(cfg);


cw=squeeze(sum(flux.data,4));
%%
figure;
imagesc(log10(abs(squeeze(cw(20,:,:)))));


title('flux, g = 0.9', 'fontsize', 20);

figure;


plot(log10(squeeze(sum(sum(cw,1),2))))


title('flux along the center, g = 0.9', 'fontsize', 20);


to see the effect, notice that I set the optical properties to match the diffusive medium benchmark used in Lihong's book as well as Fig. 6.4 in the MCML manual.

please also note that mcml allows photon trajectory to wonder far away from the launch position without truncating by a distance, while in mcx, simulation domain is bounded by your voxel size, you can minimize this boundary effect by setting a large x/y dimension sizes, but at the cost of higher memory load (not necessarily slower).

as a matter of fact, in mcx, you can directly simulate an infinite plane wave 1D decay profile in an infinite slab without needing to approximate this using radially integrated pencil beam profile. 

This is because mcx supports cyclic boundary condition. An example is provided in

https://github.com/fangq/mcx/blob/master/mcxlab/examples/demo_infinite_slab_cyclic_bc.m

I modified this script to show the similar plot to Fig. 6.4 in the MCML manual


% only clear cfg to avoid accidentally clearing other useful data


clear cfg;
cfg.nphoton=1e7;
cfg.issrcfrom0=1;
cfg.vol=uint8(ones(40,40,100));
cfg.srcdir=[0 0 1];
cfg.gpuid=1;
cfg.autopilot=1;
cfg.prop=[0 0 1 1;0.01 10 0.9 1];
cfg.tstart=0;
cfg.unitinmm=0.1;

% a uniform planar source outside the volume
cfg.srctype='planar';
cfg.srcpos=[0 0 0];
cfg.srcparam1=[40 0 0 0];
cfg.srcparam2=[0 40 0 0];
cfg.tend=5e-9;
cfg.tstep=5e-9;
cfg.bc='ccrccr';
flux=mcxlab(cfg);

fcw=flux.data*cfg.tstep;
subplot(121);
imagesc(log10(abs(squeeze(fcw(:,10,:)))))
axis equal; colorbar
set(gca,'xlim',[0 size(cfg.vol,3)]);
title('a uniform planar source incident along an infinite slab');

subplot(122);
semilogy((1:size(cfg.vol,3))*cfg.unitinmm, squeeze(mean(mean(fcw,2),1)))
title('averaged light attenuation profile');
ylabel('averaged CW fluence 1/mm^2')
xlabel('depth (mm)')


by changing the medium's n from 1 to 1.37, I was able to obtain the plots shown in the attachment.

let me know if my comments make sense to you.

Qianqian

On 10/26/23 00:46, Lukai Wang wrote:

Thank you very much Prof. Fang. I will check more carefully. Right now I think I am simulating the fluence of a pencil beam propagating in a homogeneous medium with mcml and mcx. 

My understanding of a peak of fluence at a depth of mean-free path is that photons get scattered through the propagation and the peaks appears around the depth of mean-free path, which is shown in Chapter 3 of 《Biomedical Optics: Principles and Imaging》 by Prof. Lihong Wang as follows: 

On Wednesday, October 25, 2023 at 11:05:05 PM UTC-5 Qianqian Fang wrote:

I still don't understand the rationale why you expect the fluence to have a peak at a depth of mean-free path, I also don't have an explanation why your mcml output did not show an exponential decay ... you will have to understand what you were simulating with mcml and build an equivalent one with mcx. right now, I don't think they are simulating the same thing

On 10/26/23 00:00, Lukai Wang wrote:

It is normalized fluence.

On Wednesday, October 25, 2023 at 10:58:31 PM UTC-5 Qianqian Fang wrote:

what is the y-axis of your mcml plot? the two curves have totally different shapes, making me think that you might have plotted different quantities.

On 10/25/23 23:54, Lukai Wang wrote:

Yes I used the same setting:

The domain size is 10mm x 10mm x 100mm, with voxel size of 1mm. 

The medium parameter is: μs = 0.04mm^-1, μa = 1e-3mm^-1, g = 0, n = 1.37. 

In MCX I get an exponential distribution which fits the Beer Lambert law, while I get a distribution with a maximum at mean free path in MCML:

MCX:

MCML:

Is there something wrong with my setup? Your kindly help is really appreciated.

In addition, I tried cfg.optlevel=0 in my script with my iGPUs, but the artifacts remain the same. 

Best regards,

Lukai

On Wednesday, October 25, 2023 at 12:03:25 PM UTC-5 mcx-users wrote:

another quick comment - I know many users attempted to run mcxcl/mmc with the integrated GPU because it is abundantly available - what I want to say is it is not the best hardware to take advantage of the scalability of mcxcl/mmc. even it can run correctly, its speed up over the CPU is still quite small, and mcxcl/mmc can be a lot faster (and more reliable) if you give it a discrete GPU.

if you want to help debugging this, can you set cfg.optlevel=0 in your script and see if the artifact is still there on your iGPUs? mcxcl has 4 optimization levels (see our 2018 paper), and currently the default optimization level is 1 (which enables -cl-mad-enable flag and use the native_* versions of the math functions) - I don't know if these were not implemented properly in Intel's HD/Iris runtimes, but disabling all optimization by setting cfg.optlevel=0 could tell.

Qianqian

On 10/25/23 12:49, Qianqian Fang wrote:

hi Lukai,

two comments:

1. stop using Intel HD/Iris GPUs with mcxcl - I have noticed that Intel HD/Iris GPUs had failed a subset of my built-in tests (https://github.com/fangq/mcxcl/blob/master/test/testmcx.sh), I believe these issues were caused by Intel GPU's OpenCL runtime, but haven't been able to pin-point what was the problem - all other CPUs/GPUs we have tested mcxcl have all passed, except HD/Iris integrated GPUs.

my group also purchased a discrete Intel GPU A770 earlier this year, and it also runs fine for all tests of mcxcl. so the issue seems to be specific to HD/Iris iGPUs. at this moment, please just use mcxcl on your CPU (feel free to run the testmcx.sh script in a bash shell to verify the GPU is working).

2. please do not use 2D simulations for any quantitative studies - the 2D simulation is not physically inspired - it is a faster simulation mode that I used for testing purposes (and training data generation for some of my machine-learning based works), but it is not how photons propagates in space. although MCML outputs 2D fluence output, it still simulates photon's movement in full 3D space. you should use a 3D volume with mcxcl if you want to match mcml output.

Qianqian

On 10/25/23 12:31, Lukai Wang wrote:

Dear Prof. Fang:

Thank you for the detailed response. I was using V2020 yesterday. Today I have made the following attempts with two versions of MCX (2023 and 2020) and both mcxlab (calls mcxcl.m) and mcxlabcl on my office computer (Windows) and laptop (MacOS). mcxlab and mcxlabcl always give the same result. Seems like it's openCL device that matters. Only Intel(R) Core(TM) i5-1038NG7 CPU in my Macbook seems to work. In addition, the latest version V2023 seems to give an incorrect sudden rise in flux. 

1. Windows - Intel(R) HD Graphics 530:

V.2023 / mcxlab: No artifact, but a very uncommon raise appears at around 30, which is obviously not correct. mcxlabcl gives same results.

V.2020 / mcxlab: Significant artifact. mcxlabcl gives same results.

2. MacOS - Intel(R) Iris(TM) Plus Graphics / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz

V.2023 / mcxlab / Intel(R) Iris(TM) Plus Graphics:  Uncommon raise in flux.

V.2023 / mcxlab / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz: Seems correct

V.2020 / mcxlab / Intel(R) Iris(TM) Plus Graphics: Artifact 

V.2020 / mcxlab / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz: Seems correct.

So I believe the reason for the artifact is the OpenCL device. 

Although the result by V.2020 / mcxlab / Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz is continuous and smooth. I still have some doubt, I changed anisotropy g from 0 to 0.9 in the code and the fluences seem almost the same:

But when I was comparing with Prof. Lihong Wang's code MCML (https://omlc.org/software/mc/mcml/index.html) with the same parameters, the maximum of flux appears at the depth of transport mean free path instead of the boundary and g plays an important role as follows:

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/38560103-0503-4d3e-8aac-79bf8bf20b49n%40googlegroups.com.

[Lukai Wang]

Dear Prof. Fang:

Thank you so much for the detailed response. I get your point now. 

But when I was reproduce the last code you gave me. The boundary condition does not work - the flux goes to INF. If I comment the 'cfg.bc='ccrccr'' line, the results seem right (as shown in the attached to figures). It this related to the openCL device I am using (Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz)?

Best,

Lukai

% only clear cfg to avoid accidentally clearing other useful data

clear cfg;

cfg.nphoton=1e7;

cfg.issrcfrom0=1;

cfg.vol=uint8(ones(40,40,100));

cfg.srcdir=[0 0 1];

cfg.gpuid=2;

cfg.autopilot=1;

cfg.prop=[0 0 1 1;0.01 10 0.9 1];

cfg.tstart=0;

cfg.unitinmm=0.1;

% a uniform planar source outside the volume

cfg.srctype='planar';

cfg.srcpos=[0 0 0];

cfg.srcparam1=[40 0 0 0];

cfg.srcparam2=[0 40 0 0];

cfg.tend=5e-9;

cfg.tstep=5e-9;

cfg.bc='ccrccr';

% flux=mcxlabcl(cfg);

flux=mcxlab(cfg, 'opencl');

fcw=flux.data*cfg.tstep;

subplot(121);

imagesc(log10(abs(squeeze(fcw(:,20,:)))))

axis equal; colorbar

set(gca,'xlim',[0 size(cfg.vol,3)]);

title('a uniform planar source incident along an infinite slab');

subplot(122);

semilogy((1:size(cfg.vol,3))*cfg.unitinmm, squeeze(mean(mean(fcw,2),1)))

title('averaged light attenuation profile');

ylabel('averaged CW fluence 1/mm^2')

[Qianqian Fang]

hi Lukai,

I was able to reproduce the failed cyclic bc test script on the CPU (i7-1185G7) of my laptop on both windows and Linux.

I have created an issue on github to track this problem

https://github.com/fangq/mcxcl/issues/47

although I committed a fix earlier, and it did fix the issue on the Intel CPU, I noticed that the change had broken tests on AMD GPUs, so it is not entirely fixed, and needs further debugging.

please keep an eye on this ticket. I will update once I figure out a solution.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/e093a914-3c39-4576-b7b5-cb3932e45615n%40googlegroups.com.

[Fang, Qianqian]

Hi Lukai,

 

Try download the updated mcxlabcl and mcxcl nightlybuild from https://mcx.space/nightly/win64/

 

Although mcxcl still fails the replay tests (and passed all other tests), I was able to run the previously provided scripts on windows (win11, CPU: i7-1185G7, GPU: Iris Xe), and both the Intel CPU and GPU are now giving expected results. See my attached screenshot. The updated mcxcl also passed all tests on the Intel CPU in Windows and Linux.

 

Let me know if you are able to see the same on your CPU/GPU.

 

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/d7a6316c-321f-7bb6-2f14-9495de904ef4%40neu.edu.

 

[Lukai Wang]

Dear Peof. Fang:

I downloaded the updated nightly-build version and tested on my office computer (Windows) and laptop (MacOS). Turns out it works on my laptop CPU (MacOS,  Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz):

 MacOS:

1. Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz:

2. Intel(R) Iris(TM) Plus Graphics: MATLAB breaks down every time I try.

Windows:

1. Intel(R) HD Graphics 530: 

The results seem to be correct on MacOS,  Intel(R) Core(TM) i5-1038NG7 CPU @ 2.00GHz, which is an improvement, but breakdown or seems incorrect on the other two. 

BTW, what is the replay tests you mentioned?

Best,

Lukai