support user-defined launch angle profile, fix #129

mcxlab/examples/demo_mcxlab_launchangle.m

@@ -0,0 +1,30 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
+%
+% In this example, we show how to customize phase function using cfg.invcdf
+%
+% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.sf.net
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% only clear cfg to avoid accidentally clearing other useful data
+clear cfg
+cfg.nphoton=1e7;
+cfg.vol=uint8(ones(60,60,60));
+cfg.srcpos=[30 30 1];
+cfg.srcdir=[0 0 1];
+cfg.gpuid=1;
+% cfg.gpuid='11'; % use two GPUs together
+cfg.autopilot=1;
+cfg.prop=[0 0 1 1;0.005 0.0 1 1.37];
+cfg.tstart=0;
+cfg.tend=5e-9;
+cfg.tstep=5e-9;
+cfg.isreflect=0;
+
+cfg.srctype='cone';   % source type must be a non-pencil beam, the launch angle of the selected source will be replaced by those computed by angleinvcdf
+
+% define angleinvcdf in launch angle using cfg.angleinvcdf
+cfg.angleinvcdf=linspace(1/4, 1/3);  % launch angle is uniformly distributed between [pi/4 and pi/3]
+flux=mcxlab(cfg);
+
+mcxplotvol(log10(abs(flux.data)))

src/mcx_core.h

@@ -181,6 +181,7 @@ typedef struct  __align__(16) KernelParams {
     int replaydet;                     /**< select which detector to replay, 0 for all, -1 save all separately */
     unsigned int srcnum;               /**< total number of source patterns */
     unsigned int nphase;               /**< number of samples for inverse-cdf, will be added by 2 to include -1 and 1 on the two ends */
+    unsigned int nangle;               /**< number of samples for launch angle inverse-cdf, will be added by 2 to include 0 and 1 on the two ends */
     float omega;                       /**< modulation angular frequency (2*pi*f), in rad/s, for FD/RF replay */
     unsigned char bc[12];              /**< boundary condition flags, copy the first 12 chars from cfg->bc without the terminating NULL */
 } MCXParam;

src/mcx_utils.c

@@ -330,6 +330,8 @@ void mcx_initcfg(Config* cfg) {
     cfg->gscatter = 1e9;   /** by default, honor anisotropy for all scattering, use --gscatter to reduce it */
     cfg->nphase = 0;
     cfg->invcdf = NULL;
+    cfg->nangle = 0;
+    cfg->angleinvcdf = NULL;
     memset(cfg->jsonfile, 0, MAX_PATH_LENGTH);
     memset(cfg->bc, 0, 13);
     memset(&(cfg->srcparam1), 0, sizeof(float4));
@@ -448,6 +450,10 @@ void mcx_clearcfg(Config* cfg) {
         free(cfg->invcdf);
     }
 
+    if (cfg->angleinvcdf) {
+        free(cfg->angleinvcdf);
+    }
+
     mcx_initcfg(cfg);
 }
 
@@ -2478,7 +2484,7 @@ int mcx_loadjson(cJSON* root, Config* cfg) {
     }
 
     if (Optode) {
-        cJSON* dets, *src = FIND_JSON_OBJ("Source", "Optode.Source", Optode);
+        cJSON* dets, *vv, *src = FIND_JSON_OBJ("Source", "Optode.Source", Optode);
 
         if (src) {
             subitem = FIND_JSON_OBJ("Pos", "Optode.Source.Pos", src);
@@ -2636,6 +2642,35 @@ int mcx_loadjson(cJSON* root, Config* cfg) {
                     }
                 }
             }
+
+            subitem = FIND_JSON_OBJ("AngleInverseCDF", "Optode.Source.AngleInverseCDF", src);
+
+            if (subitem) {
+                int nangle = cJSON_GetArraySize(subitem);
+                cfg->nangle = nangle + 2; /*left-/right-ends are excluded, so added 2*/
+                cfg->nangle += (cfg->nangle & 0x1); /* make cfg.nangle even number */
+
+                if (cfg->angleinvcdf) {
+                    free(cfg->angleinvcdf);
+                }
+
+                cfg->angleinvcdf = (float*)calloc(cfg->nangle, sizeof(float));
+                cfg->angleinvcdf[0] = 0.f; /*left end is always 0.f,right-end is always 1.f*/
+                vv = subitem->child;
+
+                for (i = 1; i <= nangle; i++) {
+                    cfg->angleinvcdf[i] = vv->valuedouble;
+                    vv = vv->next;
+
+                    if (cfg->angleinvcdf[i] < cfg->angleinvcdf[i - 1] || (cfg->angleinvcdf[i] > 1.f || cfg->angleinvcdf[i] < 0.f)) {
+                        MCX_ERROR(-1, "Optode.Source.AngleInverseCDF contains invalid data; it must be a monotonically increasing vector with all values between -1 and 1");
+                    }
+                }
+
+                cfg->angleinvcdf[nangle + 1] = 1.f; /*left end is always 0.f,right-end is always 1.f*/
+                cfg->angleinvcdf[cfg->nangle - 1] = 1.f;
+            }
+
         }
 
         dets = FIND_JSON_OBJ("Detector", "Optode.Detector", Optode);

src/mcx_utils.h

@@ -261,6 +261,8 @@ typedef struct MCXConfig {
     char bc[13];                 /**<boundary condition flag for [-x,-y,-z,+x,+y,+z, det(-x,-y,-z,+x,+y,+z)], last element is always NULL for string termination */
     unsigned int nphase;         /**< number of samples for inverse-cdf, will be added by 2 to include -1 and 1 on the two ends */
     float* invcdf;               /**< equal-space sampled inversion of CDF(cos(theta)) for the phase function of the zenith angle */
+    unsigned int nangle;         /**< number of samples for inverse-cdf of launch angle, will be added by 2 to include -1 and 1 on the two ends */
+    float* angleinvcdf;          /**< equal-space sampled inversion of CDF(cos(theta)) for the phase function of the zenith angle of photon launch */
 } Config;
 
 #ifdef  __cplusplus

src/mcxlab.cpp

@@ -1041,6 +1041,30 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
         cfg->invcdf[nphase + 1] = 1.f;
         cfg->invcdf[cfg->nphase - 1] = 1.f;
         printf("mcx.invcdf=[%ld];\n", cfg->nphase);
+    } else if (strcmp(name, "angleinvcdf") == 0) {
+        dimtype nangle = mxGetNumberOfElements(item);
+        double* val = mxGetPr(item);
+
+        if (cfg->angleinvcdf) {
+            free(cfg->angleinvcdf);
+        }
+
+        cfg->nangle = (unsigned int)nangle + 2;
+        cfg->nangle += (cfg->nangle & 0x1); // make cfg.nangle even number
+        cfg->angleinvcdf = (float*)calloc(cfg->nangle, sizeof(float));
+
+        for (i = 0; i < nangle; i++) {
+            cfg->angleinvcdf[i + 1] = val[i];
+
+            if (i > 0 && (val[i] < val[i - 1] || (val[i] > 1.f || val[i] < 0.f))) {
+                mexErrMsgTxt("cfg.angleinvcdf contains invalid data; it must be a monotonically increasing vector with all values between 0 and 1");
+            }
+        }
+
+        cfg->angleinvcdf[0] = 0.f;
+        cfg->angleinvcdf[nangle + 1] = 1.f;
+        cfg->angleinvcdf[cfg->nangle - 1] = 1.f;
+        printf("mcx.angleinvcdf=[%ld];\n", cfg->nangle);
     } else if (strcmp(name, "shapes") == 0) {
         int len = mxGetNumberOfElements(item);
 

src/pmcx.cpp

@@ -708,6 +708,33 @@ void parse_config(const py::dict& user_cfg, Config& mcx_config) {
         mcx_config.invcdf[mcx_config.nphase - 1] = 1.f;
     }
 
+    if (user_cfg.contains("angleinvcdf")) {
+        auto f_style_volume = py::array_t < float, py::array::f_style | py::array::forcecast >::ensure(user_cfg["angleinvcdf"]);
+
+        if (!f_style_volume) {
+            throw py::value_error("Invalid angleinvcdf field value");
+        }
+
+        auto buffer_info = f_style_volume.request();
+        unsigned int nangle = buffer_info.shape.size();
+        float* val = static_cast<float*>(buffer_info.ptr);
+        mcx_config.nangle = nangle + 2;
+        mcx_config.nangle += (mcx_config.nangle & 0x1); // make cfg.nangle even number
+        mcx_config.angleinvcdf = (float*) calloc(mcx_config.nangle, sizeof(float));
+
+        for (int i = 0; i < nangle; i++) {
+            mcx_config.angleinvcdf[i + 1] = val[i];
+
+            if (i > 0 && (val[i] < val[i - 1] || (val[i] > 1.f || val[i] < 0.f)))
+                throw py::value_error(
+                    "cfg.angleinvcdf contains invalid data; it must be a monotonically increasing vector with all values between 0 and 1");
+        }
+
+        mcx_config.angleinvcdf[0] = 0.f;
+        mcx_config.angleinvcdf[nangle + 1] = 1.f;
+        mcx_config.angleinvcdf[mcx_config.nangle - 1] = 1.f;
+    }
+
     if (user_cfg.contains("shapes")) {
         std::string shapes_string = py::str(user_cfg["shapes"]);
 

utils/mcx2json.m

@@ -60,6 +60,7 @@ function mcx2json(cfg,filestub)
 Domain=copycfg(cfg,'issrcfrom0',Domain,'OriginType',0);
 Domain=copycfg(cfg,'unitinmm',Domain,'LengthUnit');
 Domain=copycfg(cfg,'invcdf',Domain,'InverseCDF');
+Domain=copycfg(cfg,'angleinvcdf',Domain,'AngleInverseCDF');
 
 Domain.Media=cell2struct(num2cell(cfg.prop), {'mua','mus','g','n'} ,2)';