Different results from simulations with mcx and mmc

[Josef Probst]

Hi Dr. Fang,

I am currently struggling to understand why my simulation results from mcx and mmc are very different. (With a very similar simulation setup.)

I appended pictures with the simulation setup for mcx and mmc.

mmc setup:

mcx setup:

From both simulations I take the information recorded by the detector to:

- rasterize the recorded photons based on detphoton.p

- calculate the weights with: detw= mmcdetweight(detphoton,prop,1);

And from this I generated the following images:

mcx:

mmc:

In the image of the mcx simulation one can clearly see the transition between enamel and dentin whereas in the mmc simulation the center of the tooth appears the brightest.

The mcx result is very plausible, the mmc result not really.

Do you have any idea what could create such vastly different simulation outcomes?

If I should provide code or a minimal example just tell me.

Greetings

Josef Probst

[Karl-Heinz]

Hello,

I would like to add a few details:

Mr. Probst used 3D data from human teeth obtained with microcomputed tomography and segmented to separate those tissues with different optical properties (enamel, dentin, pulp). The purpose of this study is to simulate the light propagation in teeth in order to better understand diagnostic methods which use NIR light to diagnose caries. Dental enamel is considered to be rather "transparent" (little scattering, nearly no absorption), while dentin is similar to bone (scattering a lot, some absorption).

The optical properties which were used in both simulations (mcx and mmclab) are:

cfg.prop=[0 0 1 1;
          0.1 2.867 0.99 1.63;
          0.35 22.193 0.83 1.49;         
          2.8 0 1 1.333];

Sincerely

Karl-Heinz Kunzelmann

[Karl-Heinz]

Hello,

as I have no idea how to submit here example data and code I have to add a link to my github repository where I uploaded a minimal example including tooth data so that at least the mmclab part of Mr. Probst's question can be investigated.

mmclab example including geometry data

I hope this makes it easier to analyse the problem. From a dental perspective (I am a dentist) it looks like the mmclab example does not assign "enamel" properties to the outer material layer.

Sincerely

Karl-Heinz Kunzelmann

[Qianqian Fang]

from a quick look, I at least spot one issue - in your mcx simulation, the background (in the bounding box but outside of tooth) is set to 1 (with a non-air optical properties), but in the mmc mesh, the space outside of the tooth is filled with label 0.

in mcx/mmc, label 0 is treated as the "background", any photon attempts to move from a non-zero voxel/element to a zero voxel/elem will be terminated. in your mcx simulation, it looks like your photons continue to propagate/deposit energy in the bounding box, this is entirely different compared to those from mmc.

attached is the mesh plot of your mmc domain, and the corresponding optical properties are

>> cfg.prop

ans =

         0         0    1.0000    1.0000
    0.1000    2.8670    0.9900    1.6300
    0.3500   22.1930    0.8300    1.4900
    2.8000         0    1.0000    1.3330

please compare these properties and what you used for mcx, and see if there is an inconsistency/offset.

Qianqian

Karl-Heinz Kunzelmann

--
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/0e1d95b1-716b-40f3-b0e5-727b3b7fdbaen%40googlegroups.com.

[Josef Probst]

I didn't show Karl-Heinz the details of this mcx simulation, so he wasn't aware that I used a different "cfg.prop"

for mcx I used:

  cfg.prop=[0 0 1 1; % "0" voxel
            0 0 1 1; % "1" voxel
            0.1 2.867 0.99 1.63; % "2" voxel = enamel
            0.35 22.193 0.83 1.49; % "3" voxel = dentin
            2.8 0 1 1.333]; % "4" voxel = pulp

and for mmc:

  cfg.prop=[0 0 1 1; % all meshes which are generated by mmcaddsrc/det (0, -1, -2)
            0.1 2.867 0.99 1.63; % enamel
            0.35 22.193 0.83 1.49; % dentin
            2.8 0 1 1.333]; % pulp

This means that the mapping of the properties should already be correct, right? (meaning it adresses your first paragraph)

Regarding:

• in mcx/mmc, label 0 is treated as the "background", any photon attempts to move from a non-zero voxel/element to a zero voxel/elem will be terminated. in your mcx simulation, it looks like your photons continue to propagate/deposit energy in the bounding box, this is entirely different compared to those from mmc.

I am aware that this is the case for mcx. But for is this really the case for mmc? I thought that the photons get terminated as so as they leave the meshed region. And due to this one needs to mesh the region between the source and the simulation-object and between detector(s) and the simulation-object. (with mmcaddsrc/det)

And these functions create elements with "0" in the fifth column of "elem" for the space in-between and the convex-hull and "-1" or "-2" for the src/det.

I thougt that these "0"-elements behave just as defined in the first row of cfg.prop.

Also. Here is a minimal example for the mcx-simulation: http://nextcloud.mafio.si/s/idTBArRiBHjjJpk

Maybe this helps.

Thanks for your help :)

Josef

[Qianqian Fang]

On 10/1/21 6:27 AM, Josef Probst wrote:

I didn't show Karl-Heinz the details of this mcx simulation, so he wasn't aware that I used a different "cfg.prop"

for mcx I used:

  cfg.prop=[0 0 1 1; % "0" voxel
            0 0 1 1; % "1" voxel
            0.1 2.867 0.99 1.63; % "2" voxel = enamel
            0.35 22.193 0.83 1.49; % "3" voxel = dentin
            2.8 0 1 1.333]; % "4" voxel = pulp

and for mmc:

  cfg.prop=[0 0 1 1; % all meshes which are generated by mmcaddsrc/det (0, -1, -2)
            0.1 2.867 0.99 1.63; % enamel
            0.35 22.193 0.83 1.49; % dentin
            2.8 0 1 1.333]; % pulp

This means that the mapping of the properties should already be correct, right? (meaning it adresses your first paragraph)

appears so.

Regarding:

• in mcx/mmc, label 0 is treated as the "background", any photon attempts to move from a non-zero voxel/element to a zero voxel/elem will be terminated. in your mcx simulation, it looks like your photons continue to propagate/deposit energy in the bounding box, this is entirely different compared to those from mmc.

I am aware that this is the case for mcx. But for is this really the case for mmc? I thought that the photons get terminated as so as they leave the meshed region. And due to this one needs to mesh the region between the source and the simulation-object and between detector(s) and the simulation-object. (with mmcaddsrc/det)

yes, this is correct - when external detector (when label -2 appears, indicated by flag cfg->isextdet in the code) is found, we terminate the photon when it leaves the exterior mesh surface.

what I was confused was the two attachments in your original message - mmc-sim-result.png and mcx-sim-result.png, what did you plot? are these diffuse reflectance or volumetric slices? I previously thought you were plotting the volume slices, which do not make sense because I expect exact zeros in the space outside of the tooth.

also, I don't understand the big yellow circles in both images. I do not seem to see these in your mmc sample code. your mmc sample code returned a volume that have zeros in the outside as expected.

And these functions create elements with "0" in the fifth column of "elem" for the space in-between and the convex-hull and "-1" or "-2" for the src/det.

I thougt that these "0"-elements behave just as defined in the first row of cfg.prop.

Also. Here is a minimal example for the mcx-simulation: http://nextcloud.mafio.si/s/idTBArRiBHjjJpk

can't download this, because it is small, feel free to directly attach.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/db44953c-7fc7-4e98-ad80-bb576804beacn%40googlegroups.com.

[Josef Probst]

Regarding the two images. Take a look at the "%%visualisation" part of "minimal_mmc_simulation.m".

Basically take the positions where the detector was hit and create a grid based on this. Then I sort each photon into a cell of this grid and add up the weighted value.

The images are just the log() of this grid as an image.

The yellow circles are there because the source is bigger than the tooth and thus the source shines directly on the detector on some parts.

I attached the zip with my email client - no idea if this works.

[Josef Probst]

Hello,

we have now made a little progress with our problem.

A brief summary:

We have a disk light source shining through a tooth onto a wide field detector. This light source is very large at the moment, so we can see an outline of the tooth.

Our problem so far has been that the simulation with mcxlab and the simulation with mmclab do not produce matching results.

In search of an explanation, we have tested numerous hypotheses. To help you understand the problem, we have put together a minimal example and attached it as an archive. Normally we work with Octave under Linux. The programs (mmc,mcx) were compiled locally using the current files from github.

There are 6 files in the attached zip file.

1. data.mat
• volume data of the tooth
2. tooth_mcx.m
• simulation in mcx with isreflect=0 and isreflect=1
  
3. tooth_mcx_matlab_fallback.m
• same as tooth_mcx.m, but for matlab (not correct, because it does not use mcxdetweight (but good enough))
  
4. tooth_mmc.m
• simulation in mmc with isreflect=0 and isreflect=1
  
5. tooth_mmc_matlab_fallback.m
• same as tooth_mmc.m, but for matlab (not correct, because it does not use mmcdetweight(but good enough))
6. edges_from_nbins.m
• util function for visualisation in octave

So just add the paths for your environment and execute "tooth_mcx.m" and "tooth_mmc.m" (OR the matlab fallbacks if you have, no octave on hand).

Our results:

When using "cfg.isreflect = 0;" in the mmclab simulation, the detected photon distribution is very similar to that of mcxlab.
However, when using "cfg.isreflect = 1;" the result of the two simulations differs strongly.

This can be clearly seen on the attached screenshot or when running the attached simulations.

"isreflect" obviously has a large influence on the result of the mmc simulation and a small (but visible) influence in the mcx simulation.

From the documentation of isreflect, we concluded that for "isreflect=0" differences in refractive index between materials are ignored ("matched index"), while "isreflect=1" considers a refractive index mismatch.

Strangely, "isreflect=2" and "=3" produce the same images on our external detector as "=1".

Do you have any idea why isreflect can have such different effects?

Best regards to all.

Josef Probst

[Qianqian Fang]

On 10/12/21 9:04 AM, Josef Probst wrote:

Hello,

we have now made a little progress with our problem.

A brief summary:

We have a disk light source shining through a tooth onto a wide field detector. This light source is very large at the moment, so we can see an outline of the tooth.

thanks for this informative test. I am wondering if you were using your GPU to run the mmc simulation? if so, can you set cfg.gpuid=-1 and see if the CPU version outputs something more consistent?

[Josef Probst]

Thanks for the quick response.

Yes we used different GPU(s).

I just tested your suggestion, but unfortunately the results did not differ from the GPU-calculated ones.

[Qianqian Fang]

hi Josef, I was able to reproduce the images on my Linux computer (Ubuntu 20.04, Octave 5.2, latest mmc/mcx).

from the printed messages (especially absorption fraction number and detected photon number), it does look like the boundary reflection flag is effective, shown by the increased absorption from 27.6% to 42.8% - for mcx, the ratios are 27.8% and 37.25%. the reflection case has 5% difference.

one the one side, I would like you to check your detected photon accumulation code - can you take a look at the mmclab/example/demo_wide_det.m and see if you use the built-in wide-field pixel detector accumulation by setting cfg.detparam1/2 and cfg.issaveexit = 2 in your test? does it output the same as your own accumulation?

the absorption fraction difference between mcx and mmc may also give us some clue. I will look into this a bit further tomorrow - in a way, I expect the shadows of the tooth layers to be smoother when reflection is on, like in the mmc's result, because photons travel longer and get more scattered.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/bf2c1688-cc55-480d-b834-d427d508ab48n%40googlegroups.com.

[Josef Probst]

Hi Qianqian,

I checked the "accumulation code". As you can see in the attached image, the generated images look identical.

Look at the figure titles: ise=2 (issaveexit=2) uses the builtin photon accumulation. ise=1 is my implementation.

As a sidenote: I had to rotate the tooth and move src and detector, because the photon accumulation code can only handle a wide field detector in the xy-plane.

At least that's what I guessed after looking into this function: https://github.com/fangq/mmc/blob/master/src/simpmesh.c#L1081

Josef