Understanding the "replay" mode

[Zhengchen Cai]

Hi Dr.Fang,

Congratulations for the new release of the MCX and the announcing of the first MCX/MMC workshop. 

I was playing with the "replay" mode these days. It will be very convenient for us to solve the forward problem for NIRS tomography, since it provides the Jacobian matrix (sensitivity) for a pair of source and detector. To make sure we are using this function correctly, I have some questions as following (I searched in this group but there are very few questions related to this topic. They were posted at the time when "replay" model was not validated as you mentioned): 

1) Do we have a reference that gives an idea of how does "replay" model work? For the adjoint method, I've read your paper Fang Q, et al, “Combined optical Imaging and mammography of the healthy breast: optical contrast derives from breast structure and compression,”  IEEE TMI, 28(1): 30 – 42, 2009. In terms of algorithm, is replay model using this method?   

2) To get the sensitivity of a NIRS channel, we used to calculate the fluences of  the source and detector separately from mcxlab and then substitute them into the mu_a equation. In the replay mode demo, for one pair, mcx was called twice. Fluence of source was calculated first, but the detector info was also involved. In our previous implementation, we only put one light input info (either source or detector) for fluence calculation. Why do you give the detector info for the first time of calling mcx? Is it because the detp output should be specific for the detector?  

3) I don't really understand the newcfg, what are these two matrix newcfg.seed=seeds.data; newcfg.detphotons=detp.data; doing ?  Why the seed data has to be the one from the first mcx output ? Why do we need the pathlength file detp.data for replay model? 

4) What's the difference between detp.data and detp2.data? I checked the histogram of these two metric and they are quite similar. If we want to extract the pathlength info of a channel, should we use detp or detp2, I guess it should be detp2, but no idea of the reason. For me if we input source and detector both for the first and second call of mcx, the pathlength should be similar.  

5) To use the replay model for CW NIRS forward model, do we just repeat this demo for each specific channel? Each time we just update the loc and dir of source and detector the code will be same, am I right? or there are other advanced setup for replay model that is not shown in the demo?  

Thank you very much. 

Best regards,

Zhengchen

%%%%%% replay model demo%%%%%%%

clear cfg cfgs

cfg.nphoton=1e6;

cfg.vol=uint8(ones(60,60,60));

cfg.srcpos=[30 30 1];

cfg.srcdir=[0 0 1];

cfg.gpuid=1;

cfg.autopilot=1;

cfg.prop=[0 0 1 1;0.005 1 0 1.37];

cfg.tstart=0;

cfg.tend=5e-9;

cfg.tstep=5e-10;

% calculate the flux distribution with the given config

cfg.detpos=[15 30 0 2];

[flux, detp, vol, seeds]=mcxlab(cfg);   

newcfg=cfg;

newcfg.seed=seeds.data;

newcfg.outputtype='jacobian';

newcfg.detphotons=detp.data;

[flux2, detp2, vol2, seeds2]=mcxlab(newcfg);

jac=sum(flux2.data,4); 

[Qianqian Fang]

On 7/8/2017 2:09 PM, Zhengchen Cai wrote:

Hi Dr.Fang,

Congratulations for the new release of the MCX and the announcing of the first MCX/MMC workshop.

hi Zhengchen

see my reply below.

I was playing with the "replay" mode these days. It will be very convenient for us to solve the forward problem for NIRS tomography, since it provides the Jacobian matrix (sensitivity) for a pair of source and detector. To make sure we are using this function correctly, I have some questions as following (I searched in this group but there are very few questions related to this topic. They were posted at the time when "replay" model was not validated as you mentioned):

although this feature was initially implemented in mmc (and later mcx) back in 2013, we
did not have chance to thoroughly validated the method and publish our results.

https://github.com/fangq/mmc/commit/505c9e84ed9ec785076c1637250cc205ccec598e

only last year, with the help from Ruoyang Yao, a PhD student working on this project
with me and Dr. Xavier Intes at RPI, started to run some rigorous tests and validations.
We have also expanded our initial method from building mua-only Jacobians to mus.

In April, we presented this method in BODA'2017, and you can find our short conference
paper here (I attached the PDF with this email)

https://www.osapublishing.org/abstract.cfm?uri=BODA-2017-BoW3A.3

Since then, Ruoyang has been working on a journal paper, and we are about to wrap it up.
I will send you the full paper once it is published.

1) Do we have a reference that gives an idea of how does "replay" model work? For the adjoint method, I've read your paper Fang Q, et al, “Combined optical Imaging and mammography of the healthy breast: optical contrast derives from breast structure and compression,”  IEEE TMI, 28(1): 30 – 42, 2009. In terms of algorithm, is replay model using this method?  

no, replay is different from the adjoint method. it is a Monte Carlo based approach,
therefore, gives stochastic solutions. it is very fast when the detected photon is only a
small fraction of the total simulated photons; however, the down-side is the higher
stochastic noise due to small photon numbers.

in the case where the stochastic noise in the replay-derived Jacobian is not dominant,
this method only requires 1 MC background simulation from the source, and #detectors
of replays (which are typically negligible in time) to build Jacobians for all source/detector
pairs. In comparison, the adjoint method requires one to run #source + #detector
forward simulations and then multiply each pairs to form the Jacobians.

however, this is not a fair comparison because the adjoint method results in higher SNR
since both src/det solutions utilizes all simulated photons.

2) To get the sensitivity of a NIRS channel, we used to calculate the fluences of  the source and detector separately from mcxlab and then substitute them into the mu_a equation. In the replay mode demo, for one pair, mcx was called twice. Fluence of source was calculated first, but the detector info was also involved. In our previous implementation, we only put one light input info (either source or detector) for fluence calculation. Why do you give the detector info for the first time of calling mcx? Is it because the detp output should be specific for the detector? 

you should be able to find the details of the method in the attached pdf.

3) I don't really understand the newcfg, what are these two matrix newcfg.seed=seeds.data; newcfg.detphotons=detp.data; doing ?  Why the seed data has to be the one from the first mcx output ? Why do we need the pathlength file detp.data for replay model?

again, the attached paper should answer that question. To replay a photon,
we need to initialize the photon's RNG seed using the previous simulation output (seed.data).

We will also need to initialize photon's initial weight using the detp.data
last column (this is not necessary if the photons all start at weight=1. this is
only needed when using a pattered source).

4) What's the difference between detp.data and detp2.data? I checked the histogram of these two metric and they are quite similar. If we want to extract the pathlength info of a channel, should we use detp or detp2, I guess it should be detp2, but no idea of the reason. For me if we input source and detector both for the first and second call of mcx, the pathlength should be similar. 

detp.data are the detected photon info in the initial background simulation;
detp2.data are the detected photon after the replay. Because all photons should
follow their original path if their RNG seeds are set to identical values, therefore, the
two sets of detected photons should be exactly the same! we saved those
two datasets to allow us to verify that replay was successful or not.

5) To use the replay model for CW NIRS forward model, do we just repeat this demo for each specific channel? Each time we just update the loc and dir of source and detector the code will be same, am I right? or there are other advanced setup for replay model that is not shown in the demo? 

no, you can specify cfg.replaydet and set it to different detector ID,
then call mmclab for each replaydet. This will replay the photons captured
by each specified detector, thus forming the Jacobian for each pair. You can then
stack them to form the full Jacobian.

if you use mmc binary, you can use the -P/--replaydet flag to specify the
detector to be replayed. If you set it to 0, all detected photons will be replayed.
This will form a Jacobian as the summation of all src/det pairs.

let me know if you have any further questions regarding this method.

Qianqian

Thank you very much. 

Best regards,

Zhengchen

%%%%%% replay model demo%%%%%%%

clear cfg cfgs

cfg.nphoton=1e6;

cfg.vol=uint8(ones(60,60,60));

cfg.srcpos=[30 30 1];

cfg.srcdir=[0 0 1];

cfg.gpuid=1;

cfg.autopilot=1;

cfg.prop=[0 0 1 1;0.005 1 0 1.37];

cfg.tstart=0;

cfg.tend=5e-9;

cfg.tstep=5e-10;

% calculate the flux distribution with the given config

cfg.detpos=[15 30 0 2];

[flux, detp, vol, seeds]=mcxlab(cfg);   

newcfg=cfg;

newcfg.seed=seeds.data;

newcfg.outputtype='jacobian';

newcfg.detphotons=detp.data;

[flux2, detp2, vol2, seeds2]=mcxlab(newcfg);

jac=sum(flux2.data,4); 

--
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 post to this group, send email to mcx-...@googlegroups.com.
Visit this group at https://groups.google.com/group/mcx-users.
For more options, visit https://groups.google.com/d/optout.

[Zhengchen Cai]

Hi Dr.Fang,

Thank you very much for the detailed explanation and the PDF. I think I understand the general idea of "replay" mode. Just two further questions:

1) since how noisy the pMC results is rely on the number of photons that went through the "banana' shape. What would be a reasonable or minimum numbers of photons for the background simulation?  10e9 (some references mentioned this)? I noticed I only have 100 photons were detected by the detector when launching 10e5 photons. 

2) Even if adjoint method simulated fluences of detector with large number of photons, the noise level of the sensitivity of a channel should be same as replay model since SNR of Jacobians depends on the number of photons that went through the "banana' shape. In that sense, I don't understand why adjoint has higher SNR for sensitivity. I understand it has higher SNR for fluences of detector. Probably in adjoint method, the sensitivity is roughly the multiplication of two fluences, since they are both in high SNR, the production should be in high SNR. Am I right?  

One questions about the workshop, do you have any plan to record the lectures and post them in this group? It would be very nice for people who can not attend. Or at least the slides? 

Thank you so much.

Best regards,

Zhengchen

在 2017年7月8日星期六 UTC-4下午3:51:57，q.fang写道：

[Qianqian Fang]

On 07/10/2017 02:32 PM, Zhengchen Cai wrote:

Hi Dr.Fang,

Thank you very much for the detailed explanation and the PDF. I think I understand the general idea of "replay" mode. Just two further questions:

1) since how noisy the pMC results is rely on the number of photons that went through the "banana' shape. What would be a reasonable or minimum numbers of photons for the background simulation?  10e9 (some references mentioned this)? I noticed I only have 100 photons were detected by the detector when launching 10e5 photons.

the number of detected photon is linear propotional to the number of simulated ones.
if you detect 100 from 1e5 (1e5 or 10e5=1e6?) photons, then if you run 1e7 photons,
you detect something around 1e4 photons.

on the other hand, you can increase the detector radius to capture more photons.
A doubling of the radius will give you ~4x more photons.

2) Even if adjoint method simulated fluences of detector with large number of photons, the noise level of the sensitivity of a channel should be same as replay model since SNR of Jacobians depends on the number of photons that went through the "banana' shape.

that is incorrect. for adjoint, the babana shape is formed by the multiplication of
two high-SNR solutions, resulted from large photon number simulations. it is equivalent
to 2x of the total simulated photons contribute to the profile.

In that sense, I don't understand why adjoint has higher SNR for sensitivity. I understand it has higher SNR for fluences of detector. Probably in adjoint method, the sensitivity is roughly the multiplication of two fluences, since they are both in high SNR, the production should be in high SNR. Am I right? 

yes

One questions about the workshop, do you have any plan to record the lectures and post them in this group? It would be very nice for people who can not attend. Or at least the slides?

yes, all training materials, including the slides, screencasts, problem sets will
be posted on our website after the event.

Qianqian

[Zhengchen Cai]

Hi Dr.Fang,

Thank you so much. It great to know this, I wish a great success of the workshop. Have a nice day.

Best regards,

Zhengchen

在 2017年7月11日星期二 UTC-4上午10:25:40，q.fang写道：