Path Length Problem from .mch History File

417 views
Skip to first unread message

Peter McLachlan

unread,
May 5, 2015, 4:40:48 PM5/5/15
to mcx-...@googlegroups.com
Hello.

I'm looking at optical photon path lengths and I'm encountering a discrepancy.

I'll give an example. Let's say I'm simulating a semi-infinite homogeneous medium with these optical properties: mu_a = 0.01 1/mm, mu_s' = 1 1/mm, g = 0.9 and n = 1.3. I'll use a 3cm source detector separation. This is to represent the brain where I expect (from literature) to calculated a DPF (ratio of path length to source detector separation) of about 5, which corresponds to a pathlength of about 150mm.

When I run this simulation and look at the .mch history file I get a mean pathlength of more like ~425 mm. 

However, if I generate a TPSF from the flux output file and find the mean transit time, and convert it to a mean path length I get a path of about 230mm, which is much more reasonable.

Could I be doing something wrong? I can't seem to figure this out as the history file seems pretty simple in terms of its outputs but seems to give me erroneous path lengths. Is there something about the .mch file that I am missing? 

Thanks for any help.
Peter

Qianqian Fang

unread,
May 5, 2015, 4:46:33 PM5/5/15
to mcx-...@googlegroups.com
hi Peter

the correct way to calculate the mean pathlength is to use a weighted average,
weighted by each photon packet's exiting weight.

MCX does not save the packet's exiting weight, but you can quickly calculate
it by using the partial path length and mua of each medium.

did you use the above approach or just take a simple average?

Qianqian


Thanks for any help.
Peter
--
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 http://groups.google.com/group/mcx-users.
For more options, visit https://groups.google.com/d/optout.

The information in this e-mail is intended only for the person to whom it is
addressed. If you believe this e-mail was sent to you in error and the e-mail
contains patient information, please contact the Partners Compliance HelpLine at
http://www.partners.org/complianceline . If the e-mail was sent to you in error
but does not contain patient information, please contact the sender and properly
dispose of the e-mail.

Peter McLachlan

unread,
May 6, 2015, 1:36:02 PM5/6/15
to mcx-...@googlegroups.com
Hello.
I did not take a weighted average! But I was wondering if there was a weighting factor. How would I calculate the weighting factor?
Thank you so much for your reply.
Peter

Peter McLachlan

unread,
May 6, 2015, 2:51:15 PM5/6/15
to mcx-...@googlegroups.com
I'm assuming the weighting factor would just be exp(-mu_a*PL) for all media. Thank you very much.

Qianqian Fang

unread,
May 6, 2015, 3:04:30 PM5/6/15
to mcx-...@googlegroups.com
On 05/06/2015 02:51 PM, Peter McLachlan wrote:
I'm assuming the weighting factor would just be exp(-mu_a*PL) for all media. Thank you very much.

yes, that's correct. if you have more than one medium, multiply them as

exp(-mu_a_1*PL1)*exp(-mu_a_2*PL2)*...

I will add a function in the future to do this conversion, but you
should be able to write one using the above formula.

Qianqian

Peter McLachlan

unread,
May 7, 2015, 11:39:20 AM5/7/15
to mcx-...@googlegroups.com
Thank you very much.
Peter

Peter McLachlan

unread,
Jun 11, 2015, 3:46:51 PM6/11/15
to mcx-...@googlegroups.com
Hello,

I've tried taking this weighted average. Specifically, I calculate the following for each photon packet in a homogeneous medium:

PL*exp(-mu_a*PL)

And then take the mean over all the packets in the .mch history file.

This is giving me a nonsensical result, that is less than my source-detector separation. 

Am I still doing something wrong? Thanks in advance for you help.

Peter

Peter McLachlan

unread,
Jun 11, 2015, 3:51:56 PM6/11/15
to mcx-...@googlegroups.com
I seem to be off by a factor of ten. Ten times less than what I expect the path length to be.

Qianqian Fang

unread,
Jun 11, 2015, 4:17:05 PM6/11/15
to mcx-...@googlegroups.com
On 06/11/2015 03:51 PM, Peter McLachlan wrote:
I seem to be off by a factor of ten. Ten times less than what I expect the path length to be.


your weighted average seems to be incorrect.

I used mcxlab, and run the following commands:

cfg.nphoton=1e8;
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;
cfg.detpos=[30 20 1 1;30 40 1 1;20 30 1 1;40 30 1 1];   
[flux dets]=mcxlab(cfg);

% find the photons detected by det#1
idx=find(dets.data(1,:)==1);
ppl=dets.data(3,idx)';

% weighted average of pathlengths
sum(ppl.*exp(-cfg.prop(2,1).*ppl))./(sum(exp(-cfg.prop(2,1).*ppl)))

the above command gave me 79.3 mm, given the distance between
the source and detector is 10mm. This makes DPF=7.93, which looks
reasonable.

Qianqian

ALEXANDER P. Dumont

unread,
Jun 8, 2016, 10:21:58 AM6/8/16
to mcx-users
Does this give you the weight of the photons? 

Venus

unread,
Jun 16, 2016, 5:25:09 PM6/16/16
to mcx-users
So for  multi-layer medium, should it be like this? (take 2 layer)

% find the photons detected by det#1
idx=find(dets.data(1,:)==1);
ppl_1=dets.data(3,idx)'; % layer1
ppl_2=dets.data(4,idx)';  % layer2

% weighted average of pathlengths
sum((ppl_1+ppl_2).*exp(-cfg.prop(2,1).*ppl_1).*exp(-cfg.prop(3,1).*ppl_2))./(sum(exp(-cfg.prop(2,1).*ppl_1).*exp(-cfg.prop(3,1).*ppl_2)))


Is this the correct way, please?

Thank you,

Qianqian Fang

unread,
Jun 16, 2016, 11:29:42 PM6/16/16
to mcx-...@googlegroups.com
On 06/16/2016 05:25 PM, Venus wrote:
So for  multi-layer medium, should it be like this? (take 2 layer)

% find the photons detected by det#1
idx=find(dets.data(1,:)==1);
ppl_1=dets.data(3,idx)'; % layer1
ppl_2=dets.data(4,idx)';  % layer2

% weighted average of pathlengths
sum((ppl_1+ppl_2).*exp(-cfg.prop(2,1).*ppl_1).*exp(-cfg.prop(3,1).*ppl_2))./(sum(exp(-cfg.prop(2,1).*ppl_1).*exp(-cfg.prop(3,1).*ppl_2)))

this looks fine to me. although, you may simply define ppl=ppl_1+ppl_2, and use the simple formula in my previous post. it is the same when you expand it.

Qianqian

Reply all
Reply to author
Forward
0 new messages