compact implementation of ring source, close #181

fangq · Sep 3, 2023 · ad6ff4c · ad6ff4c
1 parent b0bad9f
commit ad6ff4c
Show file tree

Hide file tree

Showing 13 changed files with 87 additions and 16 deletions.
diff --git a/mcxcloud/frontend/index.html b/mcxcloud/frontend/index.html
@@ -793,7 +793,9 @@ <h3>Backend</h3>
                 "zgaussian",
                 "line",
                 "slit",
-                "pencilarray"
+                "pencilarray",
+                "hyperboloid",
+                "ring"
               ]
             },
             "Pos": {

diff --git 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
@@ -45,6 +45,8 @@ object fmSource: TfmSource
       'line'
       'slit'
       'pencilarray'
+      'hyperboloid'
+      'ring'
     )
     OnEditingDone = edSourceEditingDone
     ParentFont = False

diff --git a/pymcx b/pymcx
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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;                  /**<data zip method "zlib","gzip","base64","lzip","lzma","lz4","lz4hc"*/
     char srctype;                /**<0:pencil,1:isotropic,2:cone,3:gaussian,4:planar,5:pattern,\
-                                         6:fourier,7:arcsine,8:disk,9:fourierx,10:fourierx2d,11:zgaussian,12:line,13:slit*/
+                                         6:fourier,7:arcsine,8:disk,9:fourierx,10:fourierx2d,11:zgaussian,\
+                                         12:line,13:slit,14:pencilarray,15:pattern3d,16:hyperboloid,17:ring*/
     char outputtype;             /**<'X' output is flux, 'F' output is fluence, 'E' energy deposit*/
     char outputformat;           /**<'mc2' output is text, 'nii': binary, 'img': regular json, 'ubj': universal binary json*/
     char faststep;               /**<1 use tMCimg-like approximated photon stepping (obsolete) */

diff --git a/src/mcxlab.cpp b/src/mcxlab.cpp
@@ -892,7 +892,7 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
         int len = mxGetNumberOfElements(item);
         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", ""
                                   };
         char strtypestr[MAX_SESSION_LENGTH] = {'\0'};
 

diff --git a/src/pmcx.cpp b/src/pmcx.cpp
@@ -573,7 +573,7 @@ void parse_config(const py::dict& user_cfg, Config& mcx_config) {
         std::string src_type = py::str(user_cfg["srctype"]);
         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", ""
                                   };
         char strtypestr[MAX_SESSION_LENGTH] = {'\0'};
 

diff --git a/src/pybind11 b/src/pybind11
diff --git a/utils/mcxpreview.m b/utils/mcxpreview.m
@@ -136,12 +136,16 @@
     % rendering area-source aperature
     hsrcarea=[]; 
     if(isfield(cfg(i),'srctype'))
-        if(strcmp(cfg(i).srctype,'disk') || strcmp(cfg(i).srctype,'gaussian') || strcmp(cfg(i).srctype,'zgaussian'))
+        if(strcmp(cfg(i).srctype,'disk') || strcmp(cfg(i).srctype,'gaussian') || strcmp(cfg(i).srctype,'zgaussian') || strcmp(cfg(i).srctype,'ring'))
             if(~isfield(cfg(i),'srcparam1'))
                 error('cfg.srcparam1 is missing');
             end
             [ncyl,fcyl]=meshacylinder(srcpos,srcpos+cfg(i).srcdir(1:3)*1e-5,cfg(i).srcparam1(1)*voxelsize,0,0);
             hsrcarea=plotmesh(ncyl,fcyl{end-1},'facecolor','r','linestyle','none');
+            if(cfg(1).srcparam1(2) > 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'))