[feat] export iquv in trajectory data when cfg.istrajstokes=1

mcxlab/mcxlab.m

@@ -316,6 +316,7 @@
 %                    'M':  return photon trajectory data as the 5th output
 %                    'P':  show progress bar
 %                    'T':  save photon trajectory data only, as the 1st output, disable flux/detp/seeds outputs
+%      cfg.istrajstokes [0]: if set to 1, traj.iquv output contains the Stokes IQUV vector along trajectories
 %      cfg.maxjumpdebug: [10000000|int] when trajectory is requested in the output,
 %                     use this parameter to set the maximum position stored. By default,
 %                     only the first 1e6 positions are stored.
@@ -372,6 +373,7 @@
 %               pos: 2-4:  x/y/z/ of each trajectory position
 %                    5:    current photon packet weight
 %             srcid: 6:    source ID (>0) of the photon
+%             iquv:  7-10: Stokes IQUV vector (>0) of the photon
 %            By default, mcxlab only records the first 1e7 positions along all
 %            simulated photons; change cfg.maxjumpdebug to define a different limit.
 %
@@ -588,7 +590,10 @@
         traj.id = typecast(data(1, :), 'uint32').';
         [traj.id, idx] = sort(traj.id);
         traj.pos = traj.pos(idx, :);
-        traj.srcid = int32(data(end, :)');
+        traj.srcid = int32(data(6, :)');
+        if(size(data, 1) >= 10)
+            traj.iquv = data(7:10, :)';
+        end
         traj.data = [single(traj.id)'; data(2:end, idx)];
         newtraj(i) = traj;
     end

src/mcx_core.cu

@@ -423,17 +423,40 @@ __device__ inline void savedetphoton(float n_det[], uint* detectedphoton, float*
  * @param[in] gdebugdata: pointer to the global-memory buffer to store the trajectory info
  */
 
-__device__ inline void savedebugdata(MCXpos* p, uint id, float* gdebugdata, int srcid) {
+__device__ inline uint savedebugdata(MCXpos* p, uint id, float* gdebugdata, int srcid) {
     uint pos = atomicAdd(gjumpdebug, 1);
 
     if (pos < gcfg->maxjumpdebug) {
-        pos *= MCX_DEBUG_REC_LEN;
+        pos *= MCX_DEBUG_REC_LEN + (gcfg->istrajstokes << 2);
         ((uint*)gdebugdata)[pos++] = id;
         gdebugdata[pos++] = p->x;
         gdebugdata[pos++] = p->y;
         gdebugdata[pos++] = p->z;
         gdebugdata[pos++] = p->w;
         gdebugdata[pos++] = srcid;
+        return pos;
+    }
+
+    return 0;
+}
+
+
+/**
+ * @brief Saving photon trajectory data for debugging purposes
+ * @param[in] p: the position/weight of the current photon packet
+ * @param[in] s: the Stokes vector for polarized light simulation
+ * @param[in] id: the global index of the photon
+ * @param[in] gdebugdata: pointer to the global-memory buffer to store the trajectory info
+ */
+
+__device__ inline void savedebugstokes(MCXpos* p, Stokes* s, uint id, float* gdebugdata, int srcid) {
+    uint pos = savedebugdata(p, id, gdebugdata, srcid);
+
+    if (pos > 0 && gcfg->istrajstokes) {
+        gdebugdata[pos++] = s->i;
+        gdebugdata[pos++] = s->q;
+        gdebugdata[pos++] = s->u;
+        gdebugdata[pos++] = s->v;
     }
 }
 
@@ -1066,7 +1089,11 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*
         ppath[gcfg->partialdata] += p->w; //< sum all the remaining energy
 
         if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
-            savedebugdata(p, ((uint)f->ndone) + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]);
+            if (ispolarized && gcfg->istrajstokes) {
+                savedebugstokes(p, s, ((uint)f->ndone) + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]);
+            } else {
+                savedebugdata(p, ((uint)f->ndone) + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]);
+            }
         }
 
         if (*mediaid == 0 && *idx1d != OUTSIDE_VOLUME_MIN && *idx1d != OUTSIDE_VOLUME_MAX && gcfg->issaveref) {
@@ -1596,7 +1623,11 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*
     updateproperty<islabel, issvmc>(prop, *mediaid, t, *idx1d, media, (float3*)p, nuvox, flipdir);
 
     if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
-        savedebugdata(p, (uint)f->ndone + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[3]);
+        if (ispolarized && gcfg->istrajstokes) {
+            savedebugstokes(p, s, (uint)f->ndone + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[3]);
+        } else {
+            savedebugdata(p, (uint)f->ndone + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata, (int)ppath[3]);
+        }
     }
 
     /**
@@ -1967,7 +1998,11 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[],
                 }
 
                 if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
-                    savedebugdata(&p, (uint)f.ndone + idx * gcfg->threadphoton + umin(idx, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]);
+                    if (ispolarized && gcfg->istrajstokes) {
+                        savedebugstokes(&p, &s, (uint)f.ndone + idx * gcfg->threadphoton + umin(idx, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]);
+                    } else {
+                        savedebugdata(&p, (uint)f.ndone + idx * gcfg->threadphoton + umin(idx, gcfg->oddphotons), gdebugdata, (int)ppath[gcfg->w0offset - 1]);
+                    }
                 }
             }
 
@@ -2638,7 +2673,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) {
     size_t photoncount = 0;
 
     unsigned int printnum;
-    unsigned int tic, tic0, tic1, toc = 0, debuglen = MCX_DEBUG_REC_LEN;
+    unsigned int tic, tic0, tic1, toc = 0, debuglen = MCX_DEBUG_REC_LEN + (cfg->istrajstokes << 2);
     size_t fieldlen;
     uint3 cp0 = cfg->crop0, cp1 = cfg->crop1;
     uint2 cachebox;
@@ -2732,7 +2767,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) {
                       cp0, cp1, uint2(0, 0), cfg->minenergy, cfg->sradius* cfg->sradius, minstep* R_C0* cfg->unitinmm, cfg->srctype,
                       cfg->voidtime, cfg->maxdetphoton,
                       cfg->medianum - 1, cfg->detnum, cfg->polmedianum, cfg->maxgate, ABS(cfg->sradius + 2.f) < EPS /*isatomic*/,
-                      (uint)cfg->maxvoidstep, cfg->issaveseed > 0, (uint)cfg->issaveref, cfg->isspecular > 0,
+                      (uint)cfg->maxvoidstep, cfg->issaveseed > 0, (uint)cfg->issaveref, cfg->isspecular > 0, (uint)cfg->istrajstokes,
                       cfg->maxdetphoton * hostdetreclen, cfg->seed, (uint)cfg->outputtype, 0, 0, cfg->faststep,
                       cfg->debuglevel, cfg->savedetflag, hostdetreclen, partialdata, w0offset, cfg->mediabyte,
                       (uint)cfg->maxjumpdebug, cfg->gscatter, is2d, cfg->replaydet, cfg->srcnum,
@@ -3811,7 +3846,7 @@ is more than what your have specified (%d), please use the -H option to specify
 #ifndef MCX_CONTAINER
 
         if ((cfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) && cfg->parentid == mpStandalone && cfg->exportdebugdata) {
-            cfg->his.colcount = MCX_DEBUG_REC_LEN;
+            cfg->his.colcount = debuglen;
             cfg->his.savedphoton = cfg->debugdatalen;
             cfg->his.totalphoton = cfg->nphoton;
             cfg->his.detected = 0;

src/mcx_core.h

@@ -171,6 +171,7 @@ typedef struct  __align__(16) KernelParams {
     unsigned int issaveseed;           /**< flag if one need to save the detected photon seeds for replay */
     unsigned int issaveref;            /**< flag if one need to save diffuse reflectance data in the 0-voxel layer next to the boundary */
     unsigned int isspecular;           /**< 0 do not perform specular reflection at launch, 1 do specular reflection */
+    unsigned int istrajstokes;         /**<1 to save Stoke's vector for trajectory data only */
     unsigned int seedoffset;           /**< offset of the seed, not used */
     int seed;                          /**< RNG seed passted from the host */
     unsigned int outputtype;           /**< Type of output to be accummulated */

src/mcx_utils.c

@@ -108,7 +108,7 @@
 const char shortopt[] = {'h', 'i', 'f', 'n', 't', 'T', 's', 'a', 'g', 'b', '-', 'z', 'u', 'H', 'P',
                          'd', 'r', 'S', 'p', 'e', 'U', 'R', 'l', 'L', '-', 'I', '-', 'G', 'M', 'A', 'E', 'v', 'D',
                          'k', 'q', 'Y', 'O', 'F', '-', '-', 'x', 'X', '-', 'K', 'm', 'V', 'B', 'W', 'w', '-',
-                         '-', '-', 'Z', 'j', '-', '-', '\0'
+                         '-', '-', 'Z', 'j', '-', '-', '-', '\0'
                         };
 
 /**
@@ -128,7 +128,7 @@ const char* fullopt[] = {"--help", "--interactive", "--input", "--photon",
                          "--maxvoidstep", "--saveexit", "--saveref", "--gscatter", "--mediabyte",
                          "--momentum", "--specular", "--bc", "--workload", "--savedetflag",
                          "--internalsrc", "--bench", "--dumpjson", "--zip", "--json", "--atomic",
-                         "--srcid", ""
+                         "--srcid", "--trajstokes", ""
                         };
 
 /**
@@ -314,6 +314,7 @@ void mcx_initcfg(Config* cfg) {
     cfg->debuglevel = 0;
     cfg->issaveseed = 0;
     cfg->issaveexit = 0;
+    cfg->istrajstokes = 0;
     cfg->ismomentum = 0;
     cfg->internalsrc = 0;
     cfg->replay.seed = NULL;
@@ -1028,9 +1029,9 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c
     }
 
     if (cfg->his.detected == 0  && cfg->his.savedphoton) {
-        char colnum[] = {1, 3, 1, 1};
-        char* dtype[] = {"uint32", "single", "single", "uint32"};
-        char* dname[] = {"photonid", "p", "w0", "srcid"};
+        char colnum[] = {1, 3, 1, 1, 4};
+        char* dtype[] = {"uint32", "single", "single", "uint32", "single"};
+        char* dname[] = {"photonid", "p", "w0", "srcid", "iquv"};
         int activecol = sizeof(colnum) - (cfg->his.totalsource == 1);
         cJSON_AddItemToObject(obj, "Trajectory", dat = cJSON_CreateObject());
 
@@ -1492,6 +1493,11 @@ void mcx_preprocess(Config* cfg) {
         cfg->issavedet = 0;
     }
 
+    // if neither trajectory or polarization is enabled, disable istrajstokes flag
+    if (!(cfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) || !((cfg->mediabyte <= 4) && (cfg->polmedianum > 0))) {
+        cfg->istrajstokes = 0;
+    }
+
     for (int i = 0; i < 6; i++)
         if (cfg->bc[i] && mcx_lookupindex(cfg->bc + i, boundarycond)) {
             MCX_ERROR(-4, "unknown boundary condition specifier");
@@ -5087,6 +5093,8 @@ void mcx_parsecmd(int argc, char* argv[], Config* cfg) {
                         cfg->sradius = (isatomic) ? -2.f : 0.f;
                     } else if (strcmp(argv[i] + 2, "srcid") == 0) {
                         i = mcx_readarg(argc, argv, i, &(cfg->srcid), "int");
+                    } else if (strcmp(argv[i] + 2, "trajstokes") == 0) {
+                        i = mcx_readarg(argc, argv, i, &(cfg->istrajstokes), "int");
                     } else if (strcmp(argv[i] + 2, "internalsrc") == 0) {
                         i = mcx_readarg(argc, argv, i, &(cfg->internalsrc), "int");
                     } else {
@@ -5291,10 +5299,10 @@ int mcx_float2half2(float input[2]) {
     f2h.f[1] = input[1];
 
     /**
-    float to half conversion
-    https://stackoverflow.com/questions/3026441/float32-to-float16/5587983#5587983
-    https://gamedev.stackexchange.com/a/17410  (for denorms)
-    */
+     * float to half conversion
+     * https://stackoverflow.com/questions/3026441/float32-to-float16/5587983#5587983
+     * https://gamedev.stackexchange.com/a/17410  (for denorms)
+     */
     m = ((f2h.i[0] >> 13) & 0x03ff);
     tmp = (f2h.i[0] >> 23) & 0xff; /*exponent*/
     tmp = (tmp - 0x70) & ((unsigned int)((int)(0x70 - tmp) >> 4) >> 27);
@@ -5600,6 +5608,7 @@ where possible parameters include (the first value in [*|*] is the default)\n\
  --srcid  [0|-1,0,1,2,..]      -1 simulate multi-source separately;0 all sources\n\
                                together; a positive integer runs a single source\n\
  --internalsrc  [0|1]          set to 1 to skip entry search to speedup launch\n\
+ --trajstokes   [0|1]          set to 1 to save Stokes IQUV in trajectory data\n\
  --maxvoidstep  [1000|int]     maximum distance (in voxel unit) of a photon that\n\
                                can travel before entering the domain, if \n\
                                launched outside (i.e. a widefield source)\n\

src/mcx_utils.h

@@ -220,6 +220,7 @@ typedef struct MCXConfig {
     char issaveexit;             /**<1 save the exit position and dir of a detected photon, 0 do not save*/
     char issaveref;              /**<1 save diffuse reflectance at the boundary voxels, 0 do not save*/
     char ismomentum;             /**<1 to save momentum transfer for detected photons, implies issavedet=1*/
+    char istrajstokes;           /**<1 to save Stokes vector for trajectory data only */
     char isdumpjson;             /**<1 to save json */
     char internalsrc;            /**<1 all photons launch positions are inside non-zero voxels, 0 let mcx search entry point*/
     int  zipid;                  /**<data zip method "zlib","gzip","base64","lzip","lzma","lz4","lz4hc"*/

src/mcxlab.cpp

@@ -209,7 +209,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
      */
     for (jstruct = 0; jstruct < ncfg; jstruct++) {  /* how many configs */
         try {
-            unsigned int partialdata, hostdetreclen;
+            unsigned int partialdata, hostdetreclen, debuglen;
             printf("Running simulations for configuration #%d ...\n", jstruct + 1);
 
             /** Initialize cfg with default values first */
@@ -237,9 +237,15 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
             cfg.issavedet = (nlhs >= 2 && cfg.issavedet == 0) ? 1 : ((nlhs < 2) ? 0 : cfg.issavedet); /** save detected photon data to the 2nd output if present */
             cfg.issaveseed = (nlhs >= 4); /** save detected photon seeds to the 4th output if present */
 
+            if (nlhs >= 5 || (cfg.debuglevel & MCX_DEBUG_MOVE_ONLY)) {
+                cfg.debuglevel |= MCX_DEBUG_MOVE;
+            }
+
             /** Validate all input fields, and warn incompatible inputs */
             mcx_validatecfg(&cfg, detps, dimdetps, seedbyte);
 
+            debuglen = MCX_DEBUG_REC_LEN + (cfg.istrajstokes << 2);
+
             partialdata = (cfg.medianum - 1) * (SAVE_NSCAT(cfg.savedetflag) + SAVE_PPATH(cfg.savedetflag) + SAVE_MOM(cfg.savedetflag));
             hostdetreclen = partialdata + SAVE_DETID(cfg.savedetflag) + 3 * (SAVE_PEXIT(cfg.savedetflag) + SAVE_VEXIT(cfg.savedetflag)) + SAVE_W0(cfg.savedetflag) + 4 * SAVE_IQUV(cfg.savedetflag);
 
@@ -281,7 +287,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
             }
 
             if (nlhs >= 5 || (cfg.debuglevel & MCX_DEBUG_MOVE_ONLY)) {
-                cfg.exportdebugdata = (float*)malloc(cfg.maxjumpdebug * sizeof(float) * MCX_DEBUG_REC_LEN);
+                cfg.exportdebugdata = (float*)malloc(cfg.maxjumpdebug * sizeof(float) * debuglen);
                 cfg.debuglevel |= MCX_DEBUG_MOVE;
             }
 
@@ -324,7 +330,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
             /** if 5th output presents, output the photon trajectory data */
             if (nlhs >= 5 || (cfg.debuglevel & MCX_DEBUG_MOVE_ONLY)) {
                 int outputidx = (cfg.debuglevel & MCX_DEBUG_MOVE_ONLY) ? 0 : 4;
-                fielddim[0] = MCX_DEBUG_REC_LEN;
+                fielddim[0] = debuglen;
                 fielddim[1] = cfg.debugdatalen; // his.savedphoton is for one repetition, should correct
                 fielddim[2] = 0;
                 fielddim[3] = 0;
@@ -567,6 +573,7 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
     GET_ONE_FIELD(cfg, issaveexit)
     GET_ONE_FIELD(cfg, ismomentum)
     GET_ONE_FIELD(cfg, isspecular)
+    GET_ONE_FIELD(cfg, istrajstokes)
     GET_ONE_FIELD(cfg, replaydet)
     GET_ONE_FIELD(cfg, faststep)
     GET_ONE_FIELD(cfg, maxvoidstep)

src/pmcx.cpp

@@ -425,6 +425,7 @@ void parse_config(const py::dict& user_cfg, Config& mcx_config) {
     GET_SCALAR_FIELD(user_cfg, mcx_config, issaveexit, py::bool_);
     GET_SCALAR_FIELD(user_cfg, mcx_config, ismomentum, py::bool_);
     GET_SCALAR_FIELD(user_cfg, mcx_config, isspecular, py::bool_);
+    GET_SCALAR_FIELD(user_cfg, mcx_config, istrajstokes, py::bool_);
     GET_SCALAR_FIELD(user_cfg, mcx_config, replaydet, py::int_);
     GET_SCALAR_FIELD(user_cfg, mcx_config, faststep, py::bool_);
     GET_SCALAR_FIELD(user_cfg, mcx_config, maxvoidstep, py::int_);
@@ -1034,6 +1035,7 @@ py::dict pmcx_interface(const py::dict& user_cfg) {
     int thread_id = 0;
     size_t field_dim[6];
     py::dict output;
+    unsigned int debuglen = 0;  /** number of float numbers per photon trajectory position */
 
     try {
         /*
@@ -1043,6 +1045,8 @@ py::dict pmcx_interface(const py::dict& user_cfg) {
 
         parse_config(user_cfg, mcx_config);
 
+        debuglen = MCX_DEBUG_REC_LEN + (mcx_config.istrajstokes << 2);
+
         /** The next step, we identify gpu number and query all GPU info */
         if (!(active_dev = mcx_list_gpu(&mcx_config, &gpu_info))) {
             mcx_error(-1, "No GPU device found\n", __FILE__, __LINE__);
@@ -1094,7 +1098,7 @@ py::dict pmcx_interface(const py::dict& user_cfg) {
         }
 
         if (mcx_config.debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
-            mcx_config.exportdebugdata = (float*) malloc(mcx_config.maxjumpdebug * sizeof(float) * MCX_DEBUG_REC_LEN);
+            mcx_config.exportdebugdata = (float*) malloc(mcx_config.maxjumpdebug * sizeof(float) * debuglen);
         }
 
         /** Start multiple threads, one thread to run portion of the simulation on one CUDA GPU, all in parallel */
@@ -1131,7 +1135,7 @@ py::dict pmcx_interface(const py::dict& user_cfg) {
         field_dim[5] = 1;
 
         if (mcx_config.debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
-            field_dim[0] = MCX_DEBUG_REC_LEN;
+            field_dim[0] = debuglen;
             field_dim[1] = mcx_config.debugdatalen; // his.savedphoton is for one repetition, should correct
             field_dim[2] = 0;
             field_dim[3] = 0;