'wp' output is not consistent with detp.nscat

[andreafarina]

Dear @fangq , sorry again to bother but I'm going to checking also the Jacobian for mus now and I've found an issue. The situation is the same discussed in #91 , simple slab homogeneous. I run the reference simulation with mua = 0 to be sure that the initial weights in the replay are 1 for each photon. From the output of detp.nscat, I obtain the distribution of total scattering events by binning and histogramming the photons based on their time-of-flight.

After that, I replay the detected photons with 'wp'.

1. The first thing that is strange to me is that the data obtained are floating points. With replyweight = [1,1,..] wp, as described in your paper, should exactly be the number of scattering events in the element/node. I suppose that, due to the discretization, there is some kind of internal interpolation.
2. By summing all the wp data in the elements/nodes I rigorously have to obtain the same distribution I calculated from detp.nscat. This is not the case. wp always overestimates the total scattering events and this seems to depend both on mus and on the number of elements/nodes. It is like some scattering events are shared among the elements.

This error propagates in the final calculation of the Jacobian mus.

I attach the MATLAB code, together with the function I use for Extracting and binning. Here is the screenshot including also the comparison of wl discussed in #91 that is solved by forcing the current mua = 0 also for mua>0.

Thanks a lot for your support.

Best regards

Andrea

demo_example_replay_AF2.zip

[fangq]

hi @andreafarina, the wp output is not accumulative scattering event counts in each cell, but the weighted average - each replayed photon packet has different weight, and their detected weights have to be considered in the replay.

the raw output of wp mode is the \bar{p} item described in our replay paper, see Eq 9 in the replay paper.

to form the Jacobian of mus, one should divide \bar{p} by mus in each cell, and add the Jacobian of the mua (-\bar{L}).

we have implemented these in our provided mmcjmus.m script, please see the post-processing of the wp output here

https://github.com/fangq/mmc/blob/v2024.2/matlab/mmcjmus.m#L61-L72

let me know if this matches your expectations.

[andreafarina]

Dear Quinqian, thank you very much for the fast reply and thank you for your constant support! I’m very excited of using MMC and now I’m exploring the many possibilities of using it.

I know that, after the calculation of wp you have to post-process like in the script mmcjmus.m. I’ve also tried to compare with this and I’ve found the same discrepancies I mentioned in the GitHub issue.

Probably I forgot to mention that I use the option cfg.isnormalized = 0 like I did in ‘wl’ mode so, if the photon weights are 1 like in the case of mua = 0, in each cell I expect to find the sum of the scattering numbers of each photon passing through the cell, isn’t it correct? This is why I checked by comparing with the results of detp.nscat. For a homogeneous medium, this should be the result...

Moreover, once calculated the complete Jacobian for mus, being the slab homogeneous, summing over the cells and multiplying by a small deltaMus I expect to find the same temporal profile of TPSF(mus + dmus) - TPSF(mus). I’ve found this agreement for Jmua but not with Jmus.

Moreover, what I see is that, depending on the discretisation of the medium, this sum changes and the trend is that the discrepancy increases with the increasing of the number of elements…

By trying to debugging this on my side, I’ve performed also a couple of sanity checks:

WL newcfg.isnormalized = 0; newcfg.outputtype = 'wl'; % replay and get wl [cube2, detp2, ~, ~] = mmclab(newcfg);

newcfg.isnormalized = 1; [cube3, detp3, ~, ~] = mmclab(newcfg);

Making the elementwise ratio between cube2.data/cube3.data I obtain for each cell exactly the number of detecteod photons (not the sum of the wieghts) both with mua = 0 and mua >0. This is consistent with the bug about wl but I can easily manage it.

By performing the overall sum of wl cells and the overall sum of detp2.ppath the result must be the same, and this is the case as shown in the output:

K>> sum(cube2.data(:))

ans =

7.3136e+06

K>> sum(detp2.ppath(:))

ans =

single

 7313612

Now I do the same with wp and nscat...

WP newcfg.isnormalized = 0; newcfg.outputtype = 'wp'; % replay and get wp [cube2, detp2, ~, ~] = mmclab(newcfg);

newcfg.isnormalized = 1; [cube3, detp3, ~, ~] = mmclab(newcfg);

Making the elementwise ratio between cube2.data/cube3.data I obtain the same results of wl: the ratio is the total number of detected photons.

By performing the overall sum of wp cells and the overall sum of detp2.nscat the result must be the same, but this is not the case as shown in the output: K>> sum(detp2.nscat(:)) sum(cube2.data(:))

ans =

 7258074

ans =

 7.2988e+06

this number increases by increasing the number of elements of the mesh and it happens both with basisorder = 1 and 0. This is consistent with the issue because there is an overestimation of the scattering events in ‘wp’.

By making a very simple mesh (like you can see in the attached code) [cfg.node, cfg.elem] = genT6mesh(0:60:60, 0:60:60, 0:60:60);

I’ve observed that, with basisorder = 1, the last node value of wp is always equal to the first one, probably this can help in debugging.

I attach a revised version of the script with implemented the calculation of the Jacobian and these sanity checks…hope the comment clarify it...

I hope not to bother you too much….If you also have time or my long mail was not clear, we can organize a call with you or with a your student to clarify…

Thanks a lot again

Hope to see you in the next conference

best regards

Andrea

==========================================

Andrea Farina, PhD

Consiglio Nazionale delle Ricerche (CNR), Istituto di Fotonica e Nanotecnologie (IFN)

c/o Politecnico di Milano, Dipartimento di Fisica

Piazza Leonardo da Vinci 32

20133 Milano, Italy

Fax: +39 02 2399 6126

Tel: +39 02 2399 6024

E-mail: @. @.>

Skype account: theflours

==========================================

On 1 Sep 2024, at 18:10, Qianqian Fang @.***> wrote:

hi @andreafarina https://github.com/andreafarina, the wp output is not accumulative scattering event counts in each cell, but the weighted average - each replayed photon packet has different weight, and their detected weights have to be considered in the replay.

the raw output of wp mode is the \bar{p} item described in our replay paper, see Eq 9 in the replay paper.

to form the Jacobian of mus, one should divide \bar{p} by mus in each cell, and add the Jacobian of the mua (\bar{L}).

we have implemented these in our provided mmcjmus.m script, please see the post-processing of the wp output here

https://github.com/fangq/mmc/blob/v2024.2/matlab/mmcjmus.m#L61-L72

let me know if this matches your expectations.

— Reply to this email directly, view it on GitHub https://github.com/fangq/mmc/issues/92#issuecomment-2323410741, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADSBMMYQ3SCDIQFPNUNARB3ZUM37BAVCNFSM6AAAAABNO4JBZSVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGMRTGQYTANZUGE. You are receiving this because you were mentioned.