From ca1bf2b9da826994d7d1be12f5aa6c4ec4aa1076 Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Mon, 9 Oct 2023 01:03:21 -0400 Subject: [PATCH] use focal-length to select interpolation or discrete mode, #129 --- mcxlab/examples/demo_mcxlab_launchangle.m | 42 +++++++++++++++++++---- mcxlab/mcxlab.m | 27 +++++++++++++++ src/mcx_core.cu | 29 +++++++++++----- src/mcx_utils.c | 16 +++------ src/mcxlab.cpp | 14 +++----- src/pmcx.cpp | 15 +++----- 6 files changed, 95 insertions(+), 48 deletions(-) diff --git a/mcxlab/examples/demo_mcxlab_launchangle.m b/mcxlab/examples/demo_mcxlab_launchangle.m index da29b018..cfd79f1a 100644 --- a/mcxlab/examples/demo_mcxlab_launchangle.m +++ b/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))) \ No newline at end of file +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))) diff --git a/mcxlab/mcxlab.m b/mcxlab/mcxlab.m index 87636866..828b7444 100644 --- a/mcxlab/mcxlab.m +++ b/mcxlab/mcxlab.m @@ -112,6 +112,33 @@ % Example: % cfg.shapes: a JSON string for additional shapes in the grid % Example: +% 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 +% 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 % 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] diff --git a/src/mcx_core.cu b/src/mcx_core.cu index 5ea32000..ef9565f4 100644 --- a/src/mcx_core.cu +++ b/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; diff --git a/src/mcx_utils.c b/src/mcx_utils.c index bd9b3ece..179eed16 100644 --- a/src/mcx_utils.c +++ b/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; } } diff --git a/src/mcxlab.cpp b/src/mcxlab.cpp index a0406048..fdcf41b9 100644 --- a/src/mcxlab.cpp +++ b/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); diff --git a/src/pmcx.cpp b/src/pmcx.cpp index b1d6959e..37835cba 100644 --- a/src/pmcx.cpp +++ b/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(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(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")) {