copy detected photon post-processing scripts from mmc to mcx

utils/mcxdettime.m

@@ -0,0 +1,41 @@
+function dett=mcxdettime(detp,prop)
+%
+% dett=mcxdetweight(detp,prop)
+%
+% Recalculate the detected photon time using partial path data and 
+% optical properties (for perturbation Monte Carlo or detector readings)
+%
+% author: Qianqian Fang (q.fang <at> neu.edu)
+%	  Ruoyang Yao (yaor <at> rpi.edu) 
+%
+% input:
+%     detp: the 2nd output from mcxlab. detp can be either a struct or an array (detp.data)
+%     prop: optical property list, as defined in the cfg.prop field of mcxlab's input
+%
+% output:
+%     dett: re-caculated detected photon time 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
+%
+
+R_C0 = 3.335640951981520e-12;	% inverse of light speed in vacuum
+
+medianum=size(prop,1);
+if(medianum<=1)
+    error('empty property list');
+end
+
+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;
+    end
+else
+    detp=detp';
+    dett=zeros(size(detp,1),1);
+    for i=1:medianum-1
+        dett=dett+prop(i+1,4)*detp(:,i+medianum)*R_C0;
+    end
+end

utils/mcxdettpsf.m

@@ -0,0 +1,46 @@
+function tpsf = mcxdettpsf(detp,detnum,prop,time)
+%
+% tpsf=mcxdettpsf(detp,detnum,prop,time)
+%
+% Calculate the temporal point spread function curve of a specified detector
+% given the partial path data, optical properties, and distribution of time bins
+%
+% author: Qianqian Fang (q.fang <at> neu.edu)
+%	  Ruoyang Yao (yaor <at> rpi.edu) 
+%
+% input:
+%     detp: the 2nd output from mcxlab. detp can be either a struct or an array (detp.data)
+%     detnum: specified detector number
+%     prop: optical property list, as defined in the cfg.prop field of mcxlab's input
+%     time: distribution of time bins, a 1*3 vector [tstart tend tstep]
+%           could be defined different from in the cfg of mcxlab's input
+% output:
+%     tpsf: caculated temporal point spread function curve of the specified detector
+%
+% this file is copied from Mesh-based Monte Carlo (MMC)
+%
+% License: GPLv3, see http://mcx.space/ for details
+%
+
+% select the photon data of the specified detector
+detp.data=detp.data(:,detp.detid==detnum);
+
+% calculate the detected photon weight and arrival time
+replayweight=mcxdetweight(detp.data,prop);
+replaytime=mcxdettime(detp.data,prop);
+
+% define temporal point spread function vector
+nTG = round((time(2)-time(1))/time(3));     % maximum time gate number
+tpsf = zeros(nTG,1);
+
+% calculate the time bin, make sure not to exceed the boundary
+ntg = ceil((replaytime-time(1))/time(3));
+ntg(ntg<1)=1;
+ntg(ntg>nTG)=nTG;
+
+% add each photon weight to corresponding time bin
+for i=1:length(replayweight)
+    tpsf(ntg(i)) = tpsf(ntg(i))+replayweight(i);
+end
+
+end

utils/mcxdetweight.m

@@ -0,0 +1,56 @@
+function detw=mcxdetweight(detp,prop,w0)
+%
+% detw=mcxdetweight(detp,prop)
+%    or
+% detw=mcxdetweight(detp,prop,w0)
+%
+% Recalculate the detected photon weight using partial path data and 
+% optical properties (for perturbation Monte Carlo or detector readings)
+%
+% author: Qianqian Fang (q.fang <at> neu.edu)
+%
+% input:
+%     detp: the 2nd output from mcxlab. detp can be either a struct or an array (detp.data)
+%     prop: optical property list, as defined in the cfg.prop field of mcxlab's input
+%     w0: (optional), the initial weight of photon. 
+%           if detp is a struct, this input is ignored, detp.w0 is used instead
+%           if detp is an array and w0 is ignored, the last row of detp is used as w0
+%
+% 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
+%
+
+medianum=size(prop,1);
+if(medianum<=1)
+    error('empty property list');
+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));
+    end
+else
+    detp=detp';
+    if(nargin<3)
+        w0=detp(:,end);
+    end
+    detw=w0(:);
+    if(size(detp,2)>=2*medianum+1)
+        for i=1:medianum-1
+            detw=detw.*exp(-prop(i+1,1)*detp(:,i+medianum));
+        end
+    else
+        for i=1:medianum-1
+            detw=detw.*exp(-prop(i+1,1)*detp(:,i+2));
+        end
+    end
+end

utils/mcxmeanpath.m

@@ -0,0 +1,22 @@
+function avgpath=mcxmeanpath(detp,prop)
+%
+% avgpath=mcxmeanpath(detp,prop)
+%
+% Calculate the average pathlengths for each tissue type for a given source-detector pair
+%
+% author: Qianqian Fang (q.fang <at> neu.edu)
+%
+% input:
+%     detp: the 2nd output from mcxlab. detp can be either a struct or an array (detp.data)
+%     prop: optical property list, as defined in the cfg.prop field of mcxlab's input
+%
+% output:
+%     avepath: the average pathlength for each tissue type 
+%
+% this file is copied from Mesh-based Monte Carlo (MMC)
+%
+% License: GPLv3, see http://mcx.space/ for details
+%
+
+detw=mcxdetweight(detp,prop);
+avgpath=sum(detp.ppath.*repmat(detw(:),1,size(detp.ppath,2))) / sum(detw(:));

utils/mcxmeanscat.m

@@ -0,0 +1,22 @@
+function avgnscat=mcxmeanscat(detp,prop)
+%
+% avgnscat=mcxmeanscat(detp,prop)
+%
+% Calculate the average scattering event counts for each tissue type for a given source-detector pair
+%
+% author: Qianqian Fang (q.fang <at> neu.edu)
+%
+% input:
+%     detp: the 2nd output from mcxlab. detp can be either a struct or an array (detp.data)
+%     prop: optical property list, as defined in the cfg.prop field of mcxlab's input
+%
+% output:
+%     avgnscat: the average scattering event count for each tissue type 
+%
+% this file is copied from Mesh-based Monte Carlo (MMC)
+%
+% License: GPLv3, see http://mcx.space/ for details
+%
+
+detw=mcxdetweight(detp,prop);
+avgnscat=sum(double(detp.nscat).*repmat(detw(:),1,size(detp.nscat,2))) / sum(detw(:));