From ad6ff4cfc3c46814bdd8ebf74fb9c5fdd5ff1ab4 Mon Sep 17 00:00:00 2001 From: Qianqian Fang Date: Sun, 3 Sep 2023 10:20:21 -0400 Subject: [PATCH] compact implementation of ring source, close #181 --- mcxcloud/frontend/index.html | 4 +- mcxlab/examples/demo_mcxlab_srctype.m | 54 +++++++++++++++++++++++++++ mcxstudio/mcxsource.lfm | 2 + pymcx | 2 +- schema/mcxinput.json | 4 +- src/mcx_const.h | 1 + src/mcx_core.cu | 19 ++++++---- src/mcx_utils.c | 2 +- src/mcx_utils.h | 3 +- src/mcxlab.cpp | 2 +- src/pmcx.cpp | 2 +- src/pybind11 | 2 +- utils/mcxpreview.m | 6 ++- 13 files changed, 87 insertions(+), 16 deletions(-) diff --git a/mcxcloud/frontend/index.html b/mcxcloud/frontend/index.html index 30d1b2d7..54965e50 100644 --- a/mcxcloud/frontend/index.html +++ b/mcxcloud/frontend/index.html @@ -793,7 +793,9 @@

Backend

"zgaussian", "line", "slit", - "pencilarray" + "pencilarray", + "hyperboloid", + "ring" ] }, "Pos": { diff --git a/mcxlab/examples/demo_mcxlab_srctype.m b/mcxlab/examples/demo_mcxlab_srctype.m index 3ac4edcf..ab07e331 100644 --- a/mcxlab/examples/demo_mcxlab_srctype.m +++ b/mcxlab/examples/demo_mcxlab_srctype.m @@ -370,8 +370,62 @@ axis equal; colorbar title('a gaussian beam source'); + %% test group 5 +clear cfg; +figure; +cfg.nphoton=1e6; +cfg.vol=uint8(ones(60,60,60)); +cfg.gpuid=1; +cfg.autopilot=1; +cfg.prop=[0 0 1 1;0.005 1 0.8 1.37]; +cfg.tstart=0; +cfg.seed=99999; +cfg.tend=5e-11; +cfg.tstep=5e-11; + +% a ring beam +cfg.srctype='ring'; +cfg.srcpos=[30 30 -10]; +cfg.srcdir=[0 0 1]; +cfg.srcparam1=[15 10 0 0]; +cfg.srcparam2=[0 0 0 0]; +flux=mcxlab(cfg); +fcw=flux.data*cfg.tstep; +subplot(221); +imagesc(log10(abs(squeeze(fcw(:,:,1))))) +axis equal; colorbar +title('a ring beam'); + +% a diverging ring sector beam +cfg.srctype='ring'; +cfg.srcpos=[30 30 -10]; +cfg.srcdir=[0 0 1 -20]; +cfg.srcparam1=[15 10 0 2*pi/3]; +cfg.srcparam2=[0 0 0 0]; +flux=mcxlab(cfg); +fcw=flux.data*cfg.tstep; +subplot(222); +imagesc(log10(abs(squeeze(fcw(:,:,1))))) +axis equal; colorbar +title('a ring-sector beam'); + +% a ring beam +cfg.srctype='pencilarray'; +cfg.srcpos=[10 10 0]; +cfg.srcdir=[0 0 1]; +cfg.srcparam1=[40 0 0 4]; +cfg.srcparam2=[0 40 0 5]; +flux=mcxlab(cfg); +fcw=flux.data*cfg.tstep; +subplot(223); +imagesc(log10(abs(squeeze(fcw(:,:,1))))) +axis equal; colorbar +title('a pencil beam array)'); + +%% test group 6 + %debug flag to retrieve/test build-in RNG cfg.vol=uint8(ones(100,100,100)); cfg.debuglevel='R'; diff --git a/mcxstudio/mcxsource.lfm b/mcxstudio/mcxsource.lfm index 8a815dff..402e1b11 100644 --- a/mcxstudio/mcxsource.lfm +++ b/mcxstudio/mcxsource.lfm @@ -45,6 +45,8 @@ object fmSource: TfmSource 'line' 'slit' 'pencilarray' + 'hyperboloid' + 'ring' ) OnEditingDone = edSourceEditingDone ParentFont = False diff --git a/pymcx b/pymcx index e957869b..c5c841d0 160000 --- a/pymcx +++ b/pymcx @@ -1 +1 @@ -Subproject commit e957869b5761d76e12a8e7e665171946828c2469 +Subproject commit c5c841d02e1252c83f88c73337d6c870bb8b668c diff --git a/schema/mcxinput.json b/schema/mcxinput.json index aa34f0fd..705df7f1 100644 --- a/schema/mcxinput.json +++ b/schema/mcxinput.json @@ -216,7 +216,9 @@ "zgaussian", "line", "slit", - "pencilarray" + "pencilarray", + "hyperboloid", + "ring" ] }, "Pos": { diff --git a/src/mcx_const.h b/src/mcx_const.h index 0cd6407c..7fbdb5fb 100644 --- a/src/mcx_const.h +++ b/src/mcx_const.h @@ -93,6 +93,7 @@ #define MCX_SRC_PENCILARRAY 14 /**< a rectangular array of pencil beams */ #define MCX_SRC_PATTERN3D 15 /**< a 3D pattern source, starting from srcpos, srcparam1.{x,y,z} define the x/y/z dimensions */ #define MCX_SRC_HYPERBOLOID_GAUSSIAN 16 /**< Gaussian-beam with spot focus, scrparam1.{x,y,z} define beam waist, distance from source to focus, rayleigh range */ +#define MCX_SRC_RING 17 /**< ring/ring-sector source, scrparam1.{x,y} defines the outer/inner radius, srcparam1.{z,w} defines start/end angle*/ #define SAVE_DETID(a) ((a) & 0x1) /**< mask to save detector ID*/ #define SAVE_NSCAT(a) ((a)>>1 & 0x1) /**< output partial scattering counts */ diff --git a/src/mcx_core.cu b/src/mcx_core.cu index 3f4a023b..6fb78772 100644 --- a/src/mcx_core.cu +++ b/src/mcx_core.cu @@ -328,7 +328,7 @@ __device__ inline void updatestokes(Stokes* s, float theta, float phi, float3* u temp = (u2->z > -1.f && u2->z < 1.f) ? rsqrtf((1.f - costheta * costheta) * (1.f - u2->z * u2->z)) : 0.f; cosi = (temp == 0.f) ? 0.f : (((phi > ONE_PI && phi < TWO_PI) ? 1.f : -1.f) * (u2->z * costheta - u->z) * temp); - cosi = fmax(-1.f, fmin(cosi, 1.f)); + cosi = fmaxf(-1.f, fminf(cosi, 1.f)); sini = sqrtf(1.f - cosi * cosi); cos22 = 2.f * cosi * cosi - 1.f; @@ -1240,15 +1240,20 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime* } case (MCX_SRC_DISK): + case (MCX_SRC_RING): case (MCX_SRC_GAUSSIAN): { // uniform disk distribution or collimated Gaussian-beam // Uniform disk point picking // http://mathworld.wolfram.com/DiskPointPicking.html - float sphi, cphi; - float phi = TWO_PI * rand_uniform01(t); + float phi, sphi, cphi, r; + + if (gcfg->srctype == MCX_SRC_RING && (gcfg->srcparam1.z > 0.f || gcfg->srcparam1.w > 0.f)) { + phi = fabsf(gcfg->srcparam1.z - gcfg->srcparam1.w) * rand_uniform01(t) + fminf(gcfg->srcparam1.z, gcfg->srcparam1.w); + } else { + phi = TWO_PI * rand_uniform01(t); + } sincosf(phi, &sphi, &cphi); - float r; - if (gcfg->srctype == MCX_SRC_DISK) { + if (gcfg->srctype == MCX_SRC_DISK || gcfg->srctype == MCX_SRC_RING) { r = sqrtf(rand_uniform01(t) * fabsf(gcfg->srcparam1.x * gcfg->srcparam1.x - gcfg->srcparam1.y * gcfg->srcparam1.y) + gcfg->srcparam1.y * gcfg->srcparam1.y); } else if (fabsf(gcfg->c0.w) < 1e-5f || fabsf(gcfg->srcparam1.y) < 1e-5f) { r = sqrtf(0.5f * rand_next_scatlen(t)) * gcfg->srcparam1.x; @@ -1724,7 +1729,7 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], // in early CUDA, when ran=1, CUDA gives 1.000002 for tmp0 which produces nan later // detected by Ocelot,thanks to Greg Diamos,see http://bit.ly/cR2NMP - tmp0 = fmax(-1.f, fmin(1.f, tmp0)); + tmp0 = fmaxf(-1.f, fminf(1.f, tmp0)); theta = acosf(tmp0); stheta = sinf(theta); @@ -1854,7 +1859,7 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[], GPUDEBUG(("p=[%f %f %f] -> <%f %f %f>*%f -> hit=[%d %d %d] flip=%d\n", p.x, p.y, p.z, v.x, v.y, v.z, len, flipdir[0], flipdir[1], flipdir[2], flipdir[3])); /** if the consumed unitless scat length is less than what's left in f.pscat, keep moving; otherwise, stop in this voxel */ - slen = fmin(slen, f.pscat); + slen = fminf(slen, f.pscat); /** final length that the photon moves - either the length to move to the next voxel, or the remaining scattering length */ len = ((prop.mus == 0.f) ? len : (slen / prop.mus * (v.nscat + 1.f > gcfg->gscatter ? (1.f - prop.g) : 1.f))); diff --git a/src/mcx_utils.c b/src/mcx_utils.c index 05ad4c40..8c3aec07 100644 --- a/src/mcx_utils.c +++ b/src/mcx_utils.c @@ -196,7 +196,7 @@ const char boundarydetflag[] = {'0', '1', '\0'}; const char* srctypeid[] = {"pencil", "isotropic", "cone", "gaussian", "planar", "pattern", "fourier", "arcsine", "disk", "fourierx", "fourierx2d", "zgaussian", - "line", "slit", "pencilarray", "pattern3d", "hyperboloid", "" + "line", "slit", "pencilarray", "pattern3d", "hyperboloid", "ring", "" }; diff --git a/src/mcx_utils.h b/src/mcx_utils.h index 018424ff..ddcf44f4 100644 --- a/src/mcx_utils.h +++ b/src/mcx_utils.h @@ -211,7 +211,8 @@ typedef struct MCXConfig { char internalsrc; /**<1 all photons launch positions are inside non-zero voxels, 0 let mcx search entry point*/ int zipid; /** 0) + [ncyl,fcyl]=meshacylinder(srcpos,srcpos+cfg(i).srcdir(1:3)*1e-5,cfg(i).srcparam1(2)*voxelsize,0,0); + hsrcarea=plotmesh(ncyl,fcyl{end-1},'facecolor','k','linestyle','none'); + end elseif(strcmp(cfg(i).srctype,'planar') || strcmp(cfg(i).srctype,'pattern') || strcmp(cfg(i).srctype,'fourier') || ... strcmp(cfg(i).srctype,'fourierx') || strcmp(cfg(i).srctype,'fourierx2d') || strcmp(cfg(i).srctype,'pencilarray')) if(~isfield(cfg(i),'srcparam1') || ~isfield(cfg(i),'srcparam2'))