use focal-length to select interpolation or discrete mode, #129

mcxlab/examples/demo_mcxlab_launchangle.m

@@ -1,30 +1,58 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
 %
-% In this example, we show how to customize phase function using cfg.invcdf
+% In this example, we show how to customize photon launch angle distrbution
+% using cfg.angleinvcdf
 %
 % 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.srcpos=[30 30 15];
 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.prop=[0 0 1 1;0 0 1 1];
 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]
+cfg.angleinvcdf=linspace(1/6, 1/5, 5);  % launch angle is uniformly distributed between [pi/6 and pi/5] with interpolation (odd-number length)
 flux=mcxlab(cfg);
 
-mcxplotvol(log10(abs(flux.data)))
+mcxplotvol(log10(abs(flux.data)))
+hold on
+plot(cfg.angleinvcdf)
+
+
+%% test Lambertian/cosine distribution
+cfg.angleinvcdf=asin(0:0.01:1)/pi;  % Lambertian distribution pdf: PDF(theta)=cos(theta)/pi, CDF=sin(theta)/pi, invcdf=asin(0:1)/pi, with interpolation (cfg.srcdir(4) is not specified)
+flux=mcxlab(cfg);
+mcxplotvol(log10(abs(flux.data)))
+
+
+%% discrete angles with interpolation
+cfg.angleinvcdf=[0 0 0 0 1/6 1/6 1/4];  % interpolatation is used between discrete angular values
+flux=mcxlab(cfg);
+mcxplotvol(log10(abs(flux.data)))
+
+
+%% discrete angles without interpolation (by setting focal-length cfg.srcdir(4) to 1)
+cfg.angleinvcdf=[0 0 0 0 1/6 1/6 1/4];
+cfg.srcdir(4)=1;   % disable interpolation, use the discrete angles in angleinvcdf elements only
+flux=mcxlab(cfg);
+mcxplotvol(log10(abs(flux.data)))
+
+%% can be applied to area source - convolve with the area
+cfg.srctype='disk';
+cfg.srcparam1=[10,0,0,0];
+cfg.srcdir(4)=0;
+flux=mcxlab(cfg);
+mcxplotvol(log10(abs(flux.data)))

mcxlab/mcxlab.m

@@ -112,6 +112,33 @@
 %                      Example: <demo_sphere_cube_subpixel.m>
 %      cfg.shapes:     a JSON string for additional shapes in the grid
 %                      Example: <demo_mcxyz_skinvessel.m>
+%      cfg.invcdf:     user-specified scattering phase function. To use this, one must define
+%                      a vector with monotonically increasing value between -1 and 1
+%                      defining the discretized inverse function of the cummulative density function (CDF)
+%                      of u=cos(theta), i.e. inv(CDF(u))=inv(CDF(cos(theta))), where theta [0-pi] is the
+%                      zenith angle of the scattering event, and u=cos(theta) is between -1 and 1.
+%                      Please note that the defined inv(CDF) is relative to u, not to theta. This is because
+%                      it is relatively easy to compute as P(u) is the primary form of phase function
+%                      see <demo_mcxlab_phasefun.m>
+%      cfg.angleinvcdf: user-specified launch angle distribution. To use this, one must define
+%                      a vector with monotonically increasing value between 0 and 1
+%                      defining the discretized inverse function of the
+%                      cummulative density function (CDF) of (theta/pi),
+%                      i.e. inv(CDF(theta/pi)), where theta in the range of
+%                      [0-pi] is the zenith angle of the launch angle relative to cfg.srcdir.
+%
+%                      when source focal-length, cfg.srcdir(4), is set to 0
+%                      (default), the angleinvcdf is randomly sampled with
+%                      linear interpolation, i.e. every sample uses two
+%                      elements to interpolate the desired launch angle;
+%                      this is suited when the distribution is continuous.
+%
+%                      when cfg.srcdir(4) is set to 1, the angleinvcdf
+%                      vector is randomly sampled without interpolation
+%                      (i.e. only return the theta/pi values specified
+%                      inside this array); this is best suited for discrete
+%                      angular distribution.
+%                      see <demo_mcxlab_launchangle.m>
 %      cfg.gscatter:   after a photon completes the specified number of
 %                      scattering events, mcx then ignores anisotropy g
 %                      and only performs isotropic scattering for speed [1e9]

src/mcx_core.cu

@@ -1407,21 +1407,32 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*
                 continue;
             }
 
-            if (gcfg->nangle > 2) {
+            if (gcfg->nangle) {
                 /**
                  * If angleinvcdf is defined, use user defined launch zenith angle distribution
                  */
 
                 float ang, stheta, ctheta, sphi, cphi;
-                ang = rand_uniform01(t) * (gcfg->nangle - 2);
-                sphi = ang - ((int)ang);
-                ang = (1.f - sphi) * ((float*)(sharedmem))[((int)ang >= gcfg->nangle - 1 ? gcfg->nangle - 2 : (int)(ang)) + 1 + gcfg->nphase] +
-                      sphi * ((float*)(sharedmem))[((int)ang + 1 >= gcfg->nangle - 1 ? gcfg->nangle - 2 : (int)(ang) + 1) + 1 + gcfg->nphase];
-                ang *= ONE_PI; // next zenith angle computed based on angleinvcdf
-                sincosf(ang, &stheta, &ctheta);
+
+                if (gcfg->c0.w > 0.f) { // if focal-length > 0, no interpolation, just read the angleinvcdf value
+                    ang = fminf(rand_uniform01(t) * gcfg->nangle, gcfg->nangle - EPS);
+                    cphi = ((float*)(sharedmem))[(int)(ang) + gcfg->nphase];
+                } else { // odd number length, interpolate between neigboring values
+                    ang = fminf(rand_uniform01(t) * (gcfg->nangle - 1), gcfg->nangle - 1 - EPS);
+                    sphi = ang - ((int)ang);
+                    cphi = ((1.f - sphi) * (((float*)(sharedmem))[((int)ang >= gcfg->nangle - 1 ? gcfg->nangle - 1 : (int)(ang)) + gcfg->nphase]) +
+                            sphi * (((float*)(sharedmem))[((int)ang + 1 >= gcfg->nangle - 1 ? gcfg->nangle - 1 : (int)(ang) + 1) + gcfg->nphase]));
+                }
+
+                cphi *= ONE_PI; // next zenith angle computed based on angleinvcdf
+                sincosf(cphi, &stheta, &ctheta);
                 ang = TWO_PI * rand_uniform01(t); //next arimuth angle
                 sincosf(ang, &sphi, &cphi);
-                *((float4*)v) = gcfg->c0;
+
+                if (gcfg->c0.w < 1.5f && gcfg->c0.w >= 0.f) {
+                    *((float4*)v) = gcfg->c0;
+                }
+
                 rotatevector(v, stheta, ctheta, sphi, cphi);
             } else if (canfocus) {
                 /**
@@ -3165,7 +3176,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) {
              */
 
             /** \c ispencil: template constant, if 1, launch photon code is dramatically simplified */
-            int ispencil = (cfg->srctype == MCX_SRC_PENCIL);
+            int ispencil = (cfg->srctype == MCX_SRC_PENCIL && cfg->nangle < 2);
 
             /** \c isref: template constant, if 1, perform boundary reflection, if 0, total-absorbion boundary, can simplify kernel */
             int isref = cfg->isreflect;

src/mcx_utils.c

@@ -2335,7 +2335,6 @@ int mcx_loadjson(cJSON* root, Config* cfg) {
         if (val) {
             int nphase = cJSON_GetArraySize(val);
             cfg->nphase = nphase + 2; /*left-/right-ends are excluded, so added 2*/
-            cfg->nphase += (cfg->nphase & 0x1); /* make cfg.nphase even number */
 
             if (cfg->invcdf) {
                 free(cfg->invcdf);
@@ -2349,7 +2348,7 @@ int mcx_loadjson(cJSON* root, Config* cfg) {
                 cfg->invcdf[i] = vv->valuedouble;
                 vv = vv->next;
 
-                if (cfg->invcdf[i] < cfg->invcdf[i - 1] || (cfg->invcdf[i] > 1.f || cfg->invcdf[i] < -1.f)) {
+                if (cfg->invcdf[i] < cfg->invcdf[i - 1] || cfg->invcdf[i] > 1.f || cfg->invcdf[i] < -1.f) {
                     MCX_ERROR(-1, "Domain.InverseCDF contains invalid data; it must be a monotonically increasing vector with all values between -1 and 1");
                 }
             }
@@ -2647,28 +2646,23 @@ int mcx_loadjson(cJSON* root, Config* cfg) {
 
             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 */
+                cfg->nangle = nangle;
 
                 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++) {
+                for (i = 0; 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");
+                    if ((i > 0 && 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 0 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;
             }
 
         }

src/mcxlab.cpp

@@ -1026,19 +1026,17 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
         }
 
         cfg->nphase = (unsigned int)nphase + 2;
-        cfg->nphase += (cfg->nphase & 0x1); // make cfg.nphase even number
         cfg->invcdf = (float*)calloc(cfg->nphase, sizeof(float));
 
         for (i = 0; i < nphase; i++) {
             cfg->invcdf[i + 1] = val[i];
 
-            if (i > 0 && (val[i] < val[i - 1] || (val[i] > 1.f || val[i] < -1.f))) {
+            if ((i > 0 && val[i] < val[i - 1]) || val[i] > 1.f || val[i] < -1.f) {
                 mexErrMsgTxt("cfg.invcdf contains invalid data; it must be a monotonically increasing vector with all values between -1 and 1");
             }
         }
 
         cfg->invcdf[0] = -1.f;
-        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) {
@@ -1049,21 +1047,17 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
             free(cfg->angleinvcdf);
         }
 
-        cfg->nangle = (unsigned int)nangle + 2;
-        cfg->nangle += (cfg->nangle & 0x1); // make cfg.nangle even number
+        cfg->nangle = (unsigned int)nangle;
         cfg->angleinvcdf = (float*)calloc(cfg->nangle, sizeof(float));
 
         for (i = 0; i < nangle; i++) {
-            cfg->angleinvcdf[i + 1] = val[i];
+            cfg->angleinvcdf[i] = val[i];
 
-            if (i > 0 && (val[i] < val[i - 1] || (val[i] > 1.f || val[i] < 0.f))) {
+            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

@@ -692,19 +692,17 @@ void parse_config(const py::dict& user_cfg, Config& mcx_config) {
         unsigned int nphase = buffer_info.shape.size();
         float* val = static_cast<float*>(buffer_info.ptr);
         mcx_config.nphase = nphase + 2;
-        mcx_config.nphase += (mcx_config.nphase & 0x1); // make cfg.nphase even number
         mcx_config.invcdf = (float*) calloc(mcx_config.nphase, sizeof(float));
 
         for (int i = 0; i < nphase; i++) {
             mcx_config.invcdf[i + 1] = val[i];
 
-            if (i > 0 && (val[i] < val[i - 1] || (val[i] > 1.f || val[i] < -1.f)))
+            if ((i > 0 && val[i] < val[i - 1]) || val[i] > 1.f || val[i] < -1.f)
                 throw py::value_error(
                     "cfg.invcdf contains invalid data; it must be a monotonically increasing vector with all values between -1 and 1");
         }
 
         mcx_config.invcdf[0] = -1.f;
-        mcx_config.invcdf[nphase + 1] = 1.f;
         mcx_config.invcdf[mcx_config.nphase - 1] = 1.f;
     }
 
@@ -718,21 +716,16 @@ void parse_config(const py::dict& user_cfg, Config& mcx_config) {
         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.nangle = nangle;
         mcx_config.angleinvcdf = (float*) calloc(mcx_config.nangle, sizeof(float));
 
         for (int i = 0; i < nangle; i++) {
-            mcx_config.angleinvcdf[i + 1] = val[i];
+            mcx_config.angleinvcdf[i] = val[i];
 
-            if (i > 0 && (val[i] < val[i - 1] || (val[i] > 1.f || val[i] < 0.f)))
+            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")) {