support outputtype=length/l for saving total path lengths per voxel

see fangq/mcx@e67feed

mcxlabcl/examples/demo_mcxlab_srctype.m

@@ -49,7 +49,7 @@
 title('pencil beam launched from outside the volume');
 
 % an isotropic source outside the volume
-cfg.srctype='isotropic'
+cfg.srctype='isotropic';
 cfg.srcpos=[30 30 -10];
 cfg.tend=1e-9;
 cfg.tstep=1e-9;
@@ -190,7 +190,7 @@
 flux=mcxlabcl(cfg);
 fcw=flux.data*cfg.tstep;
 subplot(222);
-hs=slice(log10(abs(double(fcw))),[],[15 45],[1]);
+hs=slice(log10(abs(double(fcw))),[],[15 45],1);
 set(hs,'linestyle','none');
 axis equal; colorbar
 title('a slit source');

mcxlabcl/mcxlabcl.m

@@ -214,6 +214,8 @@
 % 'energy' - energy deposit per voxel
 % 'jacobian' or 'wl' - mua Jacobian (replay mode), 
 % 'nscat' or 'wp' - weighted scattering counts for computing Jacobian for mus (replay mode)
+% 'wm' - weighted momentum transfer for a source/detector pair (replay mode)
+% 'length' total pathlengths accumulated per voxel,
 % for type jacobian/wl/wp, example: <demo_mcxlab_replay.m>
 % and <demo_replay_timedomain.m>
 % cfg.session: a string for output file names (only used when no return variables)

src/mcx_core.cl

@@ -1202,9 +1202,9 @@ __kernel void mcx_main_loop(__global const uint* media,
  return;
  }
 
- float4 p = {0.f, 0.f, 0.f, -1.f}; //{x,y,z}: x,y,z coordinates,{w}:packet weight
- float4 v = gcfg->c0; //{x,y,z}: ix,iy,iz unitary direction vector, {w}:total scat event
- float4 f = {0.f, 0.f, 0.f, 0.f}; //f.w can be dropped to save register
+ float4 p = {0.f, 0.f, 0.f, -1.f}; //< {x,y,z}: x,y,z coordinates,{w}:packet weight
+ float4 v = gcfg->c0; //< {x,y,z}: ix,iy,iz unitary direction vector, {w}:total scat event
+ float4 f = {0.f, 0.f, 0.f, 0.f}; //< Photon parameter state: x/pscat: remaining scattering probability, y/t: photon elapse time, z/pathlen: total pathlen in one voxel, w/ndone: completed photons
 
  uint idx1d, idx1dold; //idx1dold is related to reflection
 
@@ -1423,6 +1423,8 @@ __kernel void mcx_main_loop(__global const uint* media,
  tshift = (int)(floor((photontof[tshift] - gcfg->twin0) * GPU_PARAM(gcfg, Rtstep))) +
  ( (GPU_PARAM(gcfg, replaydet) == -1) ? ((photondetid[tshift] - 1) * GPU_PARAM(gcfg, maxgate)) : 0);
  }
+ } else if (GPU_PARAM(gcfg, outputtype) == otL) {
+ weight = w0 * pathlen;
  }
 
  GPUDEBUG(((__constant char*)"deposit to [%d] %e, w=%f\n", idx1dold, weight, p.w));

src/mcx_host.cpp

@@ -1208,7 +1208,7 @@ is more than what your have specified (%d), please use the -H option to specify
  if (cfg->outputtype == otFluence) {
  scale[0] *= cfg->tstep;
  }
- } else if (cfg->outputtype == otEnergy) {
+ } else if (cfg->outputtype == otEnergy || cfg->outputtype == otL) {
  scale[0] = 1.f / cfg->energytot;
  } else if (cfg->outputtype == otJacobian || cfg->outputtype == otWP || cfg->outputtype == otDCS) {
  if (cfg->seed == SEED_FROM_FILE && cfg->replaydet == -1) {

src/mcx_utils.c

@@ -133,7 +133,7 @@ const char saveflag[] = {'D', 'S', 'P', 'M', 'X', 'V', 'W', '\0'};
  * p: scattering counts for computing Jacobians for mus
  */
 
-const char outputtype[] = {'x', 'f', 'e', 'j', 'p', 'm', '\0'};
+const char outputtype[] = {'x', 'f', 'e', 'j', 'p', 'm', 'l', '\0'};
 const char* outputformat[] = {"mc2", "nii", "hdr", "ubj", "tx3", "jnii", "bnii", ""};
 
 /**
@@ -4001,9 +4001,10 @@ where possible parameters include (the first value in [*|*] is the default)\n\
 \n"S_BOLD S_CYAN"\
 == Output options ==\n"S_RESET"\
  -s sessionid (--session) a string to label all output file names\n\
- -O [X|XFEJP] (--outputtype) X - output flux, F - fluence, E - energy density\n\
+ -O [X|XFEJPML](--outputtype) X - output flux, F - fluence, E - energy density\n\
  J - Jacobian (replay mode), P - scattering\n\
  event counts at each voxel (replay mode only)\n\
+ M - momentum transfer; L - total pathlength\n\
  -d [1|0] (--savedet) 1 to save photon info at detectors; 0 not save\n\
  -w [DP|DSPMXVW](--savedetflag)a string controlling detected photon data fields\n\
  /case insensitive/ 1 D output detector ID (1)\n\

src/mcx_utils.h

@@ -55,7 +55,7 @@
 #define NII_HEADER_SIZE 352
 #define GL_RGBA32F 0x8814
 
-enum TOutputType {otFlux, otFluence, otEnergy, otJacobian, otWP, otDCS}; /**< types of output */
+enum TOutputType {otFlux, otFluence, otEnergy, otJacobian, otWP, otDCS, otL}; /**< types of output */
 enum TMCXParent {mpStandalone, mpMATLAB}; /**< whether MCX is run in binary or mex mode */
 enum TOutputFormat {ofMC2, ofNifti, ofAnalyze, ofUBJSON, ofTX3, ofJNifti, ofBJNifti}; /**< output data format */
 enum TDeviceVendor {dvUnknown, dvNVIDIA, dvAMD, dvIntel, dvIntelGPU, dvAppleCPU};

src/mcxlabcl.cpp

@@ -746,7 +746,7 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
  printf("mcx.srctype='%s';\n", strtypestr);
  } else if (strcmp(name, "outputtype") == 0) {
  int len = mxGetNumberOfElements(item);
- const char* outputtype[] = {"flux", "fluence", "energy", "jacobian", "nscat", "wl", "wp", ""};
+ const char* outputtype[] = {"flux", "fluence", "energy", "jacobian", "nscat", "wl", "wp", "wm", "length", ""};
  char outputstr[MAX_SESSION_LENGTH] = {'\0'};
 
  if (!mxIsChar(item) || len == 0) {