Skip to content

Commit

Permalink
fix time/weight scaling, see mcx bug #83
Browse files Browse the repository at this point in the history
  • Loading branch information
fangq committed Nov 19, 2019
1 parent a10b743 commit a89457a
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 11 deletions.
18 changes: 12 additions & 6 deletions matlab/mmcdettime.m
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
function dett=mmcdettime(detp,prop)
function dett=mmcdettime(detp,prop,unitinmm)
%
% dett=mmcdetweight(detp,prop)
% dett=mmcdettime(detp,prop,unitinmm)
%
% Recalculate the detected photon time using partial path data and
% optical properties (for perturbation Monte Carlo or detector readings)
Expand All @@ -9,19 +9,25 @@
% Ruoyang Yao (yaor <at> rpi.edu)
%
% input:
% detp: the 2nd output from mcxlab. detp must be a struct
% detp: the 2nd output from mmclab. detp must be a struct
% prop: optical property list, as defined in the cfg.prop field of mmclab's input
% unitinmm: voxel edge-length in mm, should use cfg.unitinmm used to generate detp;
% if ignored, assume to be 1 (mm)
%
% output:
% dett: re-caculated detected photon time based on the partial path data and optical property table
%
% this file is part of Mesh-based Monte Carlo (MMC)
% this file is copied from Mesh-based Monte Carlo (MMC)
%
% License: GPLv3, see http://mcx.sf.net/mmc/ for details
% License: GPLv3, see http://mmc.space/ for details
%

R_C0 = 3.335640951981520e-12; % inverse of light speed in vacuum

if(nargin<=3)
unitinmm=1;
end

medianum=size(prop,1);
if(medianum<=1)
error('empty property list');
Expand All @@ -30,7 +36,7 @@
if(isstruct(detp))
dett=zeros(size(detp.ppath,1),1);
for i=1:medianum-1
dett=dett+prop(i+1,4)*detp.ppath(:,i)*R_C0;
dett=dett+prop(i+1,4)*detp.ppath(:,i)*R_C0*unitinmm;
end
else
error('the first input must be a struct with a subfield named "ppath"');
Expand Down
15 changes: 10 additions & 5 deletions matlab/mmcdetweight.m
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
function detw=mmcdetweight(detp,prop)
function detw=mmcdetweight(detp,prop,unitinmm)
%
% detw=mmcdetweight(detp,prop)
% detw=mmcdetweight(detp,prop,unitinmm)
%
% Recalculate the detected photon weight using partial path data and
% optical properties (for perturbation Monte Carlo or detector readings)
Expand All @@ -10,29 +10,34 @@
% input:
% detp: the 2nd output from mmclab. detp must a struct
% prop: optical property list, as defined in the cfg.prop field of mmclab's input
% unitinmm: voxel edge-length in mm, should use cfg.unitinmm used to generate detp;
% if ignored, assume to be 1 (mm)
%
% output:
% detw: re-caculated detected photon weight based on the partial path data and optical property table
%
% this file is copied from Mesh-based Monte Carlo (MMC)
%
% License: GPLv3, see http://mcx.space/ for details
% License: GPLv3, see http://mmc.space/ for details
%


medianum=size(prop,1);
if(medianum<=1)
error('empty property list');
end

if(nargin<=3)
unitinmm=1;
end

if(isstruct(detp))
if(~isfield(detp,'w0'))
detw=ones(size(detp.ppath,1),1);
else
detw=detp.w0;
end
for i=1:medianum-1
detw=detw.*exp(-prop(i+1,1)*detp.ppath(:,i));
detw=detw.*exp(-prop(i+1,1)*detp.ppath(:,i)*unitinmm);
end
else
error('the first input must be a struct with a subfield named "ppath"');
Expand Down

0 comments on commit a89457a

Please sign in to comment.