MMC - issue about "wl" output

[andreafarina]

Dear @fangq, I would like to discuss about an issue I have in checking the Jacobian mua. The situation is very simple, I started using the demo_example_replay.m. Homogeneous cube. I only add the comparison of the Jacobian with the TPSF calculated with the pathlength method with mua and mua + dmua (dmua = 1e-4), all unit are in mm and mm-1.

When the base absorption is mua = 0 and you replay the photon with replayweight = [1, 1, 1, 1....1] and isnormalized = 0 the results with "wl" seems correct. You see in each element/voxel the sum of the pathlength spent in that element. I sum up all the values in each node and the comparison with (TPSF(mua+dmua) - TPSF(mua))/dmua is correct.

When the base absorption is, e.g. mua = 0.01, the pathlengths after the reply seem amplified by a factor 1/mua, so to get the agreement with the incremental ratio you need to multiply by mua. The correct results should be pathlenghths slightly less than the one obtained in the case of mua = 0.

It seems that the starting weights of the photons are not correctly initialized when replyweight is different from [1, 1, 1, 1...1]. Can you please check this?

Thanks a lot!

Best regards

Andrea

[andreafarina]

an update: If i use the correcet weights of mua>0 but, before launching the replay with 'wl' I set the absorption property in cfg.prop(2,1) = 0, the result is correct! It seems that there is an incompatibility between the initial weight and the cfg.prop(:,1) data. Thanks again!

[andreafarina]

another issue: 'wl' gives zeros output data when cfg.method = 'grid' is set. It seems working only with "elem". Thanks a lot! Andrea

[fangq]

@andreafarina, please see the formulas I posted in https://github.com/fangq/mmc/issues/92#issuecomment-2323410741

the output wl is exactly the term \bar{L}, and wp is exactly the term \bar{p}. Their definitions are described in the screenshot in the above reply.

The formula for \bar{L} is purely related to photon's trajectories. Because mcx/mmc uses microscopic Beer-Lambert law formula (instead of those used by MCML), the trajectories of photons are only related to mus, and is independent to mua. so, I am surprised that you saw a dependency to mua. Could it be possible that you used MCML's postprocessing to handle this output?

also, please take a look at our mmcjmua.m script, this should directly give you the Jacobian of mua, without needing additional processing

https://github.com/fangq/mmc/blob/v2024.2/matlab/mmcjmua.m

let me know if this gives you what you wanted.

[andreafarina]

Dear Qinqian, this is exactly what I expect…the trajectories are independent on mus and I would expect not to see a dependence on mua… I tried to flag cfg.mcmethod = 0 and 1 but things remain pretty the same. I post here a test code to show what I obtain. I'm replaying forcing the replayweight = [1,1,1...1] but I change the absorption of the medium starting from 0 to 1mm-1. Then I extract an element of the output after reply and you can see the dependence on the absorption. In principle, these numbers should be exactly the same...

clear cfg newcfg;

cfg.nphoton = 1e5;
cfg.seed = 1648335518;

[cfg.node, cfg.elem] = genT6mesh(0:60:60, 0:60:60, 0:60:60);
cfg.elemprop = ones(size(cfg.elem, 1), 1);
cfg.srcpos = [30.1, 30.2, 0];
cfg.srcdir = [0 0 1];
cfg.tstart = 0;
cfg.tend = 2e-9;
cfg.tstep = 5e-11;
cfg.prop = [0 0 1 1; 
            0.00 0.10 0.01 1.0];

cfg.debuglevel = 'TP';%TP';
cfg.method = 'elem';%'elem';
%cfg.gpuid = 1;
cfg.basisorder = 0;
cfg.isreflect = 0;
cfg.detpos = [30. 20. 0. 1.];%

cfg.issaveexit = 1;  % save detected photon exit position and angles
cfg.issaveseed = 1;  % save detected photon seeds to replay later
cfg.mcmethod = 0;  % as in mmclab.m 0:MMC 1:MCML

mua = [0.0, 0.01, 0.05 1];
for i = 1:numel(mua)
    newcfg = mmclab(cfg, 'prep');  % preprocessing of the mesh to get the missing fields
    newcfg.prop(2,1) = mua(i);
    [cube, detp, ncfg, seeds] = mmclab(newcfg); % initial simulation

    newcfg.seed = seeds.data;
    % calculate the detected photon weight using the partial path output and prop
    newcfg.replayweight = mmcdetweight(detp, newcfg.prop);

    % force weight to be 1
    newcfg.replayweight = ones(size(newcfg.replayweight));

    newcfg.replaytime = mmcdettime(detp, newcfg.prop);

    newcfg.isnormalized = 0;
    newcfg.outputtype = 'wl';   % replay and get wl

    [cube2, detp2] = mmclab(newcfg);
    % the two detected photon arrays should be the same. however, because
    % the program uses multi-threading, the orders may be different
    if (all(ismember(round(detp.data' * 1e10) * 1e-10, round(detp2.data' * 1e10) * 1e-10, 'rows')))
        disp('replay is successful :-)');
    else
        disp('replay failed :-(');
    end
    one_elem_wl(i) = cube2.data(1,3); % see an element of the output
end
disp(num2str(one_elem_wl))

Output is:

312.53772 30516.6953 7005.77393 571.356201

Concerning the mmcjmua.m script after I managed a bug in detp.data handling for many sources, I've obtained exactly the same results of my processing. When you set an absorption > 0 the values of the Jacobian, normalized or unnormalized, are rougly amplified of a factor 1/mua...

If you look at the code I lastly appended in #92 you can test both wl and wp with the perturbed TPSF. Thank you very much again Best regards

Andrea

[andreafarina]

I've found the problematic line of code. https://github.com/fangq/mmc/blob/c67142db035e3a58bf0a7b460dd03bcb5421a7fa/src/mmc_core.cl#L705C17-L705C91 Adding the test condition on WL the problem is solved. I managed to compile the code under a linux workstation and it works.

if (prop.mua > 0.f) {
                if (GPU_PARAM(gcfg, outputtype) != otEnergy && GPU_PARAM(gcfg, outputtype) != otWP && GPU_PARAM(gcfg, outputtype) != otWL) {
                    ww /= prop.mua;
                }
            }