Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Estimation of the pathlengh #132

Closed
Edouard2laire opened this issue Dec 8, 2021 · 2 comments
Closed

Estimation of the pathlengh #132

Edouard2laire opened this issue Dec 8, 2021 · 2 comments
Labels

Comments

@Edouard2laire
Copy link

Edouard2laire commented Dec 8, 2021

Hello,

I am not sure if this is the correct place to ask questions. I just wanted to check if what I was doing with mcxlab was correct. Here I am trying to use it to estimate the pathlength that should be used in the modified beer lambert law.

My questions are:

  • To increase the reliability of the estimation of the mean pathlength is it better to increase the number of photons used, or increase the number of respin (or is it equivalent)? Does the duration of the simulation (tstart, tstep,stend affect this estimation?)

  • if I look at the distribution of pathlength in detphoton.ppath i get the following:
    pathlengh_distribution
    Is it normal to have photons traveling in the white matter? The estimate of the mean pathlengh isn't really aligned with the peak of the distribution, is that normal? i guess it is related to the weighting applied when computing the mean or because the distribution isn't a gaussian distribution but maybe more a log-normal one?

  • If I want to compare the pathlengh estimated by mcxlab to a fixed value, I guess it could be easier to have access to not only the mean pathlengh but the full distribution, is it possible? ( i guess detphoton.ppath isn't directly what I am interested in since we need to apply some weighting ie mcxdetweight ). Or would you have some suggestions to see if the estimated value is different from a fixed value?

  • is it possible to visualize the 'banana' shape going from the source to the detector? similarly to this image
    one_layer_expl

  • Is it possible to use multiple channels in one simulation 9eg 1 source and 2 detectors for example?) When I tried it, as soon as I have more than one detector, I get an error saying out of memory.

Sorry for the very long message :)

PS:
I am using the following configuration to start the simulation?

mcx.gpuid=1;
mcx.autopilot=1;
mcx.respin=1;
mcx.seed=1648335518;
mcx.nphoton=1e+08;
mcx.dim=[256 256 192];
mcx.mediabyte=4;
mcx.unitinmm=1;
mcx.isreflect=1;
mcx.isrefint=1;
mcx.tstart=0;
mcx.tend=5e-09;
mcx.tstep=5e-09;
mcx.isnormalized=1;
mcx.issrcfrom0=0;
mcx.medianum=6;
mcx.srcpos=[178 192 104];
mcx.srcdir=[-0.782397 -0.590269 -0.198591 0];
mcx.minenergy=0;
mcx.detnum=1;
@fangq fangq added the question label Dec 8, 2021
@fangq
Copy link
Owner

fangq commented Dec 8, 2021

To increase the reliability of the estimation of the mean pathlength is it better to increase the number of photons used, or increase the number of respin (or is it equivalent)? Does the duration of the simulation (tstart, tstep,stend affect this estimation?)

you can use one of 3 ways to run a simulation with a desired photon number N:

  1. a single mcx/mcxlab call with nphoton=N, or
  2. a single mcx session but abs(cfg.respin) number of kernel calls, with each kernel simulating N/respin photons, or
  3. running mcx for M times, with each one running N/M photons

the above approaches are ranked in the order of most robust/fastest/most convenient to the least order.

also, approach 1 has the least overhead - in approach 2, kernels have to be called multiple times, adding overhead of host-gpu data exchange; for approach 3, additional overhead such as file IO are added. Approach 2 is supposed to automatically taking care of the data normalization, while in 3, you will have to take care of averaging by yourself.

I added 2 as a workaround to the "kernel launch timed out error" that commonly seen in mcx, so a single kernel execution won't exceed the watchdog time limit. Mode 3 is largely used in distributed systems.

if your GPU is set up to run more a longer simulation, I would avoid using 2/3. The only downside is that when you are running a single session with very high photon number, say 10^10, then, there is a risk of floating-point overflow (see #41) near your source region, if you observe that, you should split into smaller sessions, which avoids the overflow problem.

Is it normal to have photons traveling in the white matter? The estimate of the mean pathlengh isn't really aligned with the peak of the distribution, is that normal? i guess it is related to the weighting applied when computing the mean or because the distribution isn't a gaussian distribution but maybe more a log-normal one?

yes, it is normal. Each partial path follows different distributions - I had previously hypothesized Erlang distribution of different parameters, but I have never had chance to prove it. It makes the most sense to me because summation of multiple exponential distribution (which the unitless scattering length follows) follows Erlang distrbution:

https://en.wikipedia.org/wiki/Erlang_distribution#Related_distributions

If I want to compare the pathlengh estimated by mcxlab to a fixed value, I guess it could be easier to have access to not only the mean pathlengh but the full distribution, is it possible? ( i guess detphoton.ppath isn't directly what I am interested in since we need to apply some weighting ie mcxdetweight ). Or would you have some suggestions to see if the estimated value is different from a fixed value?

yes, people use total pathlengths to estimate DPF (diferential pathlength factor), see https://www.nature.com/articles/pr19962544

partial path is tricky because it depends on the geometry.

is it possible to visualize the 'banana' shape going from the source to the detector? similarly to this image
one_layer_expl

yes, you can use either replay or adjoint, see replay sample script here
https://github.com/fangq/mcx/blob/master/mcxlab/examples/demo_mcxlab_replay.m

the adjoint is even simpler - you just run forward at your source, and also at your detector, then multiply the two solutions, see our replay paper Eq 14
https://www.osapublishing.org/boe/fulltext.cfm?uri=boe-9-10-4588&id=396723

Is it possible to use multiple channels in one simulation 9eg 1 source and 2 detectors for example?) When I tried it, as soon as I have more than one detector, I get an error saying out of memory.

I would like to see the error message (and your GPU model),
we never had issues

https://github.com/fangq/mcx/blob/master/mcxlab/examples/demo_colin27_atlas.m#L40-L43

@fangq
Copy link
Owner

fangq commented Dec 13, 2021

@Edouard2laire, let me know if I can close this ticket or you have additional questions.

@fangq fangq closed this as completed Jan 25, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

2 participants