Issues with 2-layer skin simulation (epidermis/dermis) + detector output questions

[Jonas]

Hi,

I am currently trying to simulate a simplified skin model consisting of two tissue layers (epidermis and dermis).
My goal is to measure the reflected photons using a detector in order to estimate the "color" of the tissue.
For this purpose, I run 3 separate simulations for three wavelengths (R, G, B), and for each wavelength I simulate two different “"lood levels" by modifying the optical properties of the dermis.

However, I am running into two issues:

1. Photons do not seem to reach the second layer (dermis) 

Although the geometry appears correct, the detector shows no photons in ppath entering the second layer.
I have verified the mesh and layer definitions, but I may still be missing something.

2. Detector output structure  [flux, det] = mmclab(cfg); 
I am unsure how to interpret the fields inside det.data.
What exactly does det.data contain (detector_id, weight?, path lengths in the tissue?)?

Thanks in advance! 

This is my code:

x_size = 500;

y_size = 500;

z_epidermis = 10;

z_dermis = 1000;

cfg.unitinmm = 0.01; % 1 grid unit equals 0.01mm (10mu m)

cfg.nphoton = 1e5; % n photons

cfg.srctype = 'disk';

cfg.srcparam1 = [100 0 0 0]; % 2 grid unit radius

cfg.srcpos = [x_size/2 y_size/2 0]; % Source on surface

cfg.srcdir = [0 0 1];

cfg.detpos = [x_size/2, y_size/2, 0, 10000];

cfg.tstart = 0;

cfg.tend = 5e-9;

cfg.tstep = 5e-9;

cfg.issavepr = 1;

cfg.issodet = 1;

cfg.outputtype = 'energy';

if not(load_file)

% Define Mesh

% Layer 1: Epidermis (0.1mm thick -> 10 units)

% Layer 2: Dermis (Infinite -> big...)

[n, f] = latticegrid([0 x_size], [0 y_size], [0 z_epidermis z_dermis]);

seeds = [x_size/2 y_size/2 z_epidermis/2; x_size/2 y_size/2 z_dermis/2]; % Seed 1: Epidermis, Seed 2: Dermis

[cfg.node, cfg.elem] = s2m(n, f, 1, 100, 'tetgen', seeds, []);

figure(1);

plotmesh(cfg.node, cfg.elem)

ticks = xticks;

xticklabels(string(ticks * cfg.unitinmm));

yticklabels(string(yticks * cfg.unitinmm));

zticklabels(string(zticks * cfg.unitinmm));

xlabel('x (mm)');

ylabel('y (mm)');

zlabel('z (mm)');

cfg.elemprop = cfg.elem(:,5);

cfg.elem = cfg.elem(:,1:4);

end

save("simulation.mat", "cfg");

%% Optical Properties

% 650nm, 550nm, 450nm

wavelengths = {"Red", "Green", "Blue"};

% blood fractions in percent

blood_fractions = [0.01, 0.05];

% melanin fractions in percent

f_melanin = 0.06;

% Epidermis [mua, mus, g, n] in 1/mm

props.epidermis.Red = [1.59, 9.41, 0.8, 1.40];

props.epidermis.Green = [2.86, 4.722, 0.85, 1.40];

props.epidermis.Blue = [4.76, 2.598, 0.9, 1.40];

% Bloodless Dermis

props.dermis.Red = [0.01, 38.69, 0.85, 1.37];

props.dermis.Green = [0.18, 48.01, 0.85, 1.37];

props.dermis.Blue = [0.1, 62.22, 0.85, 1.37];

% Skinbaseline

props.skinbaseline.Red = [0.0295, 9.41, 0.85, 1.37];

props.skinbaseline.Green = [0.0459, 4.722, 0.85, 1.37];

props.skinbaseline.Blue = [0.1219, 2.598, 0.85, 1.37];

% Oxygenated Blood

props.blood.Red = [0.2165, 10];

props.blood.Green = [25.3059, 10];

props.blood.Blue = [36.9529, 10];

props.melanin.Red = 28.3493;

props.melanin.Green = 49.4466;

props.melanin.Blue = 96.4598;

%% Simulation Loop

% define how many loops per wavelength to simulate

loops = [1, 1, 1]; % Photons (loops * cfg.nphotons)

detect_weight = zeros(2, 3);

energy_starting = zeros(2, 3);

for s = 1:2 % loop over blood fraction states

f_blood = blood_fractions(s);

for w = 1:3 % loop over wavelengths

wv = wavelengths{w};

absorbed = [];

launched_energy = [];

total_detected_energy = 0;

mua_e = f_melanin * props.melanin.(wv)(1) + (1-f_melanin) * props.skinbaseline.(wv)(1);

mua_d = f_blood * props.blood.(wv)(1) + (1-f_blood) * props.skinbaseline.(wv)(1)

mus_d = f_blood * props.blood.(wv)(2) + (1-f_blood) * props.skinbaseline.(wv)(2);

cfg.prop = [

0 0 1 1; % ID 0

mua_e, props.epidermis.(wv)(2:3), 1.0; % ID 1

mua_d, mus_d, props.dermis.(wv)(3), props.dermis.(wv)(4) % ID 2

];

total_detected_weight = 0;

total_launched = 0;

for i = 1:loops(w) % loops

[flux, det] = mmclab(cfg);

batch_weights = mmcdetweight(det, cfg.prop, cfg.unitinmm);

total_detected_weight = total_detected_weight + sum(batch_weights);

total_launched = total_launched + cfg.nphoton;

end

detect_weight(s, w) = total_detected_weight;

energy_starting(s, w) = total_launched;

end

end

save("result.mat", "det");

disp("Finished simulation!")

[Qianqian Fang]

hi Jonas,

first, your code could not be copied and executed, because of unexpected line wraps. I had to make multiple edits in order to run some of it.

Going over your code, I found two problems, one is related to mmc, one is related to your cfg setting.

first, I noticed that mmc currently disables cfg.unitinmm if running on the GPU. Without the scaling of unitinmm, the epidermis layer is over 10 mm long and it is extremely rare for photon to pass to the 2nd layer.

more specifically, it disables cfg.unitinmm when using DMMC (dual-grid MMC, which output the fluence as a 3D array like mcx), which is the default on the GPU. I can't remember why this was disabled, likely because it could accidentally cause excessive memory allocation and cause crash.

for this issue, I created a new issue on github (https://github.com/fangq/mmc/issues/117) and was able to fix it with this patch

https://github.com/fangq/mmc/commit/2a58c43376d5c6700a56a33ace67ca852831b30f

secondly, I noticed a bigger problem is that you used a distributed/widefield source (srctype='disk') without applying mesh re-tessellation. It crashes matlab.

a disk source, like planar/pattern, launches photons from distributed locations. when a photon is launched at a different location, the index of the enclosing element is different. Handling such source is significantly more complex compared to a point source (pencil, isotropic, cone). We had to do an extra step, called mesh re-tessellation to modify the mesh in order to ensure the entire source area is enclosed inside a single tetrahedron. The details of this method is described in our paper

https://www.osapublishing.org/boe/abstract.cfm?uri=boe-7-1-171

see a similar discussion on this:

https://groups.google.com/g/mmc-users/c/__1qBrhYsAg/m/RXsZk6zBAQAJ

Qianqian

On 12/10/25 09:44, Jonas wrote:

You don't often get email from jonas.k...@gmail.com. Learn why this is important

--
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/c2a53a9c-7724-44ca-804e-e70b9e102de7n%40googlegroups.com.

[Qianqian Fang]

2. Detector output structure  [flux, det] = mmclab(cfg); 
I am unsure how to interpret the fields inside det.data.
What exactly does det.data contain (detector_id, weight?, path lengths in the tissue?)?

just to answer your 2nd question: 

yes, det.data is the concatenation of the other fields in the det struct, please see the format in this section of the help info

https://github.com/fangq/mmc/blob/v2025.10/mmclab/mmclab.m#L201-L213

however, you are not recommended to use det.data. 

It is the raw output, and is sbsequently split into individual det subfields cast to the proper data types. you should use the split data struct.

On 12/10/25 09:44, Jonas wrote:

You don't often get email from jonas.k...@gmail.com. Learn why this is important

[Jonas]

Thank you for your answers.

I am trying to compile MMC from the source you updated on windows so I can use it in MATLAB. I followed the instructions in the README and created a mexops.bat pointing to Cygwin and x86_64-w64-mingw32-g++. However, when i run make mex in the src folder i only get errors:

Supported compiler not detected. You can install the freely available MinGW-w64
C/C++ compiler; visit https://www.mathworks.com/matlabcentral/fileexchange/52848
-matlab-support-for-mingw-w64-c-c-compiler. For more options, visit https://www.
mathworks.com/support/compilers.
make: *** [../commons/Makefile_common.mk:294: ../bin/mmc] Error 127

I don't know if you can help me with that...