support --bc or cfg.bc to set an entire bounding face as detector

mcxlab/examples/demo_bc_det.m

@@ -0,0 +1,29 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang
+%
+% In this example, we show how to set the 7-12 elements of cfg.bc to
+% indicate a planar detector along a boundary facet
+%
+% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.sf.net
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clear cfg cfgs
+cfg.nphoton=1e5;
+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 1 0 1.37];
+cfg.tstart=0;
+cfg.tend=5e-9;
+cfg.tstep=1e-10;
+
+% calculate the flux distribution with the given config
+cfg.bc='______111110';  % capture photons existing from all faces except z=z_max
+cfg.savedetflag='dpx';
+
+[flux,detpt]=mcxlab(cfg);
+plot3(detpt.p(:,1),detpt.p(:,2),detpt.p(:,3),'r.');
+view([0 0 -1])

mcxlab/mcxlab.m

@@ -83,6 +83,14 @@
 %                      'a': like cfg.isreflect=0, total absorption BC
 %                      'm': mirror or total reflection BC
 %                      'c': cyclic BC, enter from opposite face
+%
+%                      in addition, cfg.bc can contain up to 12 characters,
+%                      with the 7-12 characters indicating bounding box
+%                      facets -x,-y,-z,+x,+y,+z are used as a detector. The 
+%                      acceptable characters for digits 7-12 include
+%                      '0': this face is not used to detector photons
+%                      '1': this face is used to capture photons (if output detphoton)
+%                      see <demo_bc_det.m>
 %      cfg.isnormalized:[1]-normalize the output fluence to unitary source, 0-no reflection
 %      cfg.isspecular: 1-calculate specular reflection if source is outside, [0] no specular reflection
 %      cfg.maxgate:    the num of time-gates per simulation

src/mcx_const.h

@@ -38,6 +38,7 @@
 #define NO_LAUNCH          9999                    /**< when fail to launch, for debug */
 #define OUTSIDE_VOLUME_MIN 0xFFFFFFFF              /**< flag indicating the index is outside of the volume from x=xmax,y=ymax,z=zmax*/
 #define OUTSIDE_VOLUME_MAX 0x7FFFFFFF              /**< flag indicating the index is outside of the volume from x=0/y=0/z=0*/
+#define BOUNDARY_DET_MASK  0xFFFF0000              /**< flag indicating a boundary face is used as a detector*/
 #define MAX_PROP_AND_DETECTORS   4000              /**< maximum number of property + number of detectors */
 #define SEED_FROM_FILE      -999                   /**< special flag indicating to read seeds from an mch file for replay */
 

src/mcx_core.cu

@@ -181,8 +181,9 @@ __device__ inline void saveexitppath(float n_det[],float *ppath,MCXpos *p0,uint
  * @param[in] seeddata: the RNG seed of the photon at launch, need to save for replay 
  */
 
-__device__ inline void savedetphoton(float n_det[],uint *detectedphoton,float *ppath,MCXpos *p0,MCXdir *v,RandType t[RAND_BUF_LEN],RandType *seeddata,uint *idx1d){
-      uint detid=finddetector(p0);
+__device__ inline void savedetphoton(float n_det[],uint *detectedphoton,float *ppath,MCXpos *p0,MCXdir *v,RandType t[RAND_BUF_LEN],RandType *seeddata,uint isdet){
+      int detid;
+      detid=(isdet==OUTSIDE_VOLUME_MIN)?-1:(int)finddetector(p0);
       if(detid){
 	 uint baseaddr=atomicAdd(detectedphoton,1);
 	 if(baseaddr<gcfg->maxdetphoton){
@@ -766,7 +767,7 @@ __device__ inline int launchnewphoton(MCXpos *p,MCXdir *v,MCXtime *f,float3* rv,
       // let's handle detectors here
           if(gcfg->savedet){
              if((isdet&DET_MASK)==DET_MASK && *mediaid==0 && gcfg->issaveref<2)
-	         savedetphoton(n_det,dpnum,ppath,p,v,photonseed,seeddata,idx1d);
+	         savedetphoton(n_det,dpnum,ppath,p,v,photonseed,seeddata,isdet);
              clearpath(ppath,gcfg->partialdata);
           }
 #endif
@@ -1335,9 +1336,9 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
           if(p.x<0.f||p.y<0.f||p.z<0.f||p.x>=gcfg->maxidx.x||p.y>=gcfg->maxidx.y||p.z>=gcfg->maxidx.z){
               /** if photon moves outside of the volume, set mediaid to 0 */
 	      mediaid=0;
-	      isdet=gcfg->bc[(!(p.x<0.f||p.y<0.f||p.z<0.f))*3+flipdir];  /** isdet now stores the boundary condition flag, this will be overwriten before the end of the loop */
-              GPUDEBUG(("moving outside: [%f %f %f], idx1d [%d]->[out], bcflag %d\n",p.x,p.y,p.z,idx1d,isdet));
 	      idx1d=(p.x<0.f||p.y<0.f||p.z<0.f) ? OUTSIDE_VOLUME_MIN : OUTSIDE_VOLUME_MAX;
+	      isdet=gcfg->bc[(idx1d==OUTSIDE_VOLUME_MAX)*3+flipdir];  /** isdet now stores the boundary condition flag, this will be overwriten before the end of the loop */
+              GPUDEBUG(("moving outside: [%f %f %f], idx1d [%d]->[out], bcflag %d\n",p.x,p.y,p.z,idx1d,isdet));
 	  }else{
               /** otherwise, read the optical property index */
 	      mediaid=media[idx1d];
@@ -1439,8 +1440,9 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
 		 }
 	      }
               GPUDEBUG(("direct relaunch at idx=[%d] mediaid=[%d], ref=[%d] bcflag=%d timegate=%d\n",idx1d,mediaid,gcfg->doreflect,isdet,f.t>gcfg->twin1));
-	      if(launchnewphoton<ispencil, isreflect, islabel>(&p,&v,&f,&rv,&prop,&idx1d,field,&mediaid,&w0,(mediaidold & DET_MASK),ppath,
-	          n_det,detectedphoton,t,(RandType*)(sharedmem+threadIdx.x*gcfg->issaveseed*RAND_BUF_LEN*sizeof(RandType)),
+	      if(launchnewphoton<ispencil, isreflect, islabel>(&p,&v,&f,&rv,&prop,&idx1d,field,&mediaid,&w0,
+	          (((idx1d==OUTSIDE_VOLUME_MAX && gcfg->bc[9+flipdir]) || (idx1d==OUTSIDE_VOLUME_MIN && gcfg->bc[6+flipdir]))? OUTSIDE_VOLUME_MIN : (mediaidold & DET_MASK)),
+	          ppath, n_det,detectedphoton,t,(RandType*)(sharedmem+threadIdx.x*gcfg->issaveseed*RAND_BUF_LEN*sizeof(RandType)),
 		  media,srcpattern,idx,(RandType*)n_seed,seeddata,gdebugdata,gprogress))
                    break;
               isdet=mediaid & DET_MASK;
@@ -1494,7 +1496,8 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
                         transmit(&v,n1,prop.n,flipdir);
                         if(mediaid==0){ // transmission to external boundary
                             GPUDEBUG(("transmit to air, relaunch\n"));
-		    	    if(launchnewphoton<ispencil, isreflect, islabel>(&p,&v,&f,&rv,&prop,&idx1d,field,&mediaid,&w0,(mediaidold & DET_MASK),
+		    	    if(launchnewphoton<ispencil, isreflect, islabel>(&p,&v,&f,&rv,&prop,&idx1d,field,&mediaid,&w0,
+			        (((idx1d==OUTSIDE_VOLUME_MAX && gcfg->bc[9+flipdir]) || (idx1d==OUTSIDE_VOLUME_MIN && gcfg->bc[6+flipdir]))? OUTSIDE_VOLUME_MIN : (mediaidold & DET_MASK)),
 			        ppath,n_det,detectedphoton,t,(RandType*)(sharedmem+threadIdx.x*gcfg->issaveseed*RAND_BUF_LEN*sizeof(RandType)),
 				media,srcpattern,idx,(RandType*)n_seed,seeddata,gdebugdata,gprogress))
                                 break;
@@ -2021,7 +2024,7 @@ void mcx_run_simulation(Config *cfg,GPUInfo *gpu){
          param.idx1dorig=(int(floorf(p0.z))*dimlen.y+int(floorf(p0.y))*dimlen.x+int(floorf(p0.x)));
          param.mediaidorig=(cfg->vol[param.idx1dorig] & MED_MASK);
      }
-     memcpy(&(param.bc),cfg->bc,8);
+     memcpy(&(param.bc),cfg->bc,12);
      Vvox=cfg->steps.x*cfg->steps.y*cfg->steps.z; /*Vvox: voxel volume in mm^3*/
      if(cfg->seed>0)
      	srand(cfg->seed+threadid);

src/mcx_core.h

@@ -148,7 +148,7 @@ typedef struct  __align__(16) KernelParams {
   unsigned int is2d;                 /**< is the domain a 2D slice? */
   int replaydet;                     /**< select which detector to replay, 0 for all, -1 save all separately */
   unsigned int srcnum;               /**< total number of source patterns */
-  unsigned char bc[8];               /**< boundary conditions */
+  unsigned char bc[12];              /**< boundary conditions */
 }MCXParam;
 
 void mcx_run_simulation(Config *cfg,GPUInfo *gpu);

src/mcx_utils.c

@@ -166,6 +166,14 @@ const char *outputformat[]={"mc2","nii","hdr","ubj","tx3","jnii","bnii",""};
 
 const char boundarycond[]={'_','r','a','m','c','\0'};
 
+/**
+ * Boundary detection flags
+ * 0: do not detect photon
+ * 1: detect photon at that boundary
+ */
+
+const char boundarydetflag[]={'0','1','\0'};
+
 /**
  * Source type specifier
  * User can specify the source type using a string
@@ -287,7 +295,7 @@ void mcx_initcfg(Config *cfg){
      cfg->dx=cfg->dy=cfg->dz=NULL;
      cfg->gscatter=1e9;     /** by default, honor anisotropy for all scattering, use --gscatter to reduce it */
      memset(cfg->jsonfile,0,MAX_PATH_LENGTH);
-     memset(cfg->bc,0,8);
+     memset(cfg->bc,0,12);
      memset(&(cfg->srcparam1),0,sizeof(float4));
      memset(&(cfg->srcparam2),0,sizeof(float4));
      memset(cfg->deviceid,0,MAX_DEVICE);
@@ -1142,6 +1150,19 @@ void mcx_writeconfig(char *fname, Config *cfg){
  */
 
 void mcx_prepdomain(char *filename, Config *cfg){
+     int isbcdet=0;
+
+     for(int i=0;i<6;i++)
+        if(cfg->bc[i]>='A' && mcx_lookupindex(cfg->bc+i,boundarycond))
+	   MCX_ERROR(-4,"unknown boundary condition specifier");
+
+     for(int i=6;i<12;i++){
+        if(cfg->bc[i]>='0' && mcx_lookupindex(cfg->bc+i,boundarydetflag))
+	   MCX_ERROR(-4,"unknown boundary detection flags");
+	if(cfg->bc[i])
+	   isbcdet=1;
+     }
+
      if(cfg->isdumpjson==2){
 	  mcx_savejdata(cfg->jsonfile, cfg);
 	  exit(0);
@@ -1164,7 +1185,7 @@ void mcx_prepdomain(char *filename, Config *cfg){
 		mcx_convertrow2col(&(cfg->vol), &(cfg->dim));
 		cfg->isrowmajor=0;
 	}
-        if(cfg->issavedet && cfg->detnum==0)
+        if(cfg->issavedet && cfg->detnum==0 && isbcdet)
             cfg->issavedet=0;
         if(cfg->issavedet==0){
             cfg->issaveexit=0;
@@ -1205,10 +1226,6 @@ void mcx_prepdomain(char *filename, Config *cfg){
         if(cfg->deviceid[i]=='0')
            cfg->deviceid[i]='\0';
 
-     for(int i=0;i<6;i++)
-        if(cfg->bc[i]>='A' && mcx_lookupindex(cfg->bc+i,boundarycond))
-	   MCX_ERROR(-4,"unknown boundary condition specifier");
-
      if((cfg->mediabyte==MEDIA_AS_F2H || cfg->mediabyte==MEDIA_MUA_FLOAT || cfg->mediabyte==MEDIA_AS_HALF) && cfg->medianum<2)
          MCX_ERROR(-4,"the 'prop' field must contain at least 2 rows for the requested media format");
      if((cfg->mediabyte==MEDIA_ASGN_BYTE || cfg->mediabyte==MEDIA_AS_SHORT) && cfg->medianum<3)
@@ -2773,7 +2790,7 @@ void mcx_parsecmd(int argc, char* argv[], Config *cfg){
 		     	        break;
                      case 'B':
                                 if(i<argc+1)
-				    strncpy(cfg->bc,argv[i+1],8);
+				    strncpy(cfg->bc,argv[i+1],12);
 				i++;
                                	break;
 		     case 'd':
@@ -3174,6 +3191,14 @@ where possible parameters include (the first value in [*|*] is the default)\n\
 			       'a': like -b 0, total absorption BC\n\
 			       'm': mirror or total reflection BC\n\
 			       'c': cyclic BC, enter from opposite face\n\
+\n\
+			       if input contains additional 6 letters,\n\
+			       the 7th-12th letters can be:\n\
+			       '0': do not use this face to detect photon, or\n\
+			       '1': use this face for photon detection (-d 1)\n\
+			       the order of the faces for letters 7-12 is \n\
+			       the same as the first 6 letters\n\
+			       eg: --bc ______010 saves photons exiting at y=0\n\
  -u [1.|float] (--unitinmm)    defines the length unit for the grid edge\n\
  -U [1|0]      (--normalize)   1 to normalize flux to unitary; 0 save raw\n\
  -E [0|int|mch](--seed)        set random-number-generator seed, -1 to generate\n\

src/mcx_utils.h

@@ -227,7 +227,7 @@ typedef struct MCXConfig{
 	float *exportdebugdata;      /**<pointer to the buffer where the photon trajectory data are stored*/
 	uint mediabyte;              /**< how many bytes per media index, mcx supports 1, 2 and 4, 4 is the default*/
 	float *dx, *dy, *dz;         /**< anisotropic voxel spacing for x,y,z axis */
-	char bc[8];                  /**<boundary condition flag for [-x,-y,-z,+x,+y,+z,unused,unused] */
+	char bc[12];                 /**<boundary condition flag for [-x,-y,-z,+x,+y,+z, det(-x,-y,-z,+x,+y,+z)] */
 } Config;
 
 #ifdef	__cplusplus

src/mcxlab.cpp

@@ -712,7 +712,7 @@ void mcx_set_field(const mxArray *root,const mxArray *item,int idx, Config *cfg)
         jsonshapes[len]='\0';
     }else if(strcmp(name,"bc")==0){
         int len=mxGetNumberOfElements(item);
-        if(!mxIsChar(item) || len==0 || len>7)
+        if(!mxIsChar(item) || len==0 || len>12)
              mexErrMsgTxt("the 'bc' field must be a non-empty string");
 
         mxGetString(item, cfg->bc, len+1);
@@ -855,8 +855,9 @@ void mcx_replay_prep(Config *cfg){
  */
 
 void mcx_validate_config(Config *cfg){
-     int i,gates,idx1d;
+     int i,gates,idx1d, isbcdet=0;
      const char boundarycond[]={'_','r','a','m','c','\0'};
+     const char boundarydetflag[]={'0','1','\0'};
      unsigned int partialdata=(cfg->medianum-1)*(SAVE_NSCAT(cfg->savedetflag)+SAVE_PPATH(cfg->savedetflag)+SAVE_MOM(cfg->savedetflag));
      unsigned int hostdetreclen=partialdata+SAVE_DETID(cfg->savedetflag)+3*(SAVE_PEXIT(cfg->savedetflag)+SAVE_VEXIT(cfg->savedetflag))+SAVE_W0(cfg->savedetflag);
 
@@ -912,6 +913,13 @@ void mcx_validate_config(Config *cfg){
         if(cfg->bc[i]>='A' && mcx_lookupindex(cfg->bc+i,boundarycond))
 	   mexErrMsgTxt("unknown boundary condition specifier");
 
+     for(i=6;i<12;i++){
+        if(cfg->bc[i]>='0' && mcx_lookupindex(cfg->bc+i,boundarydetflag))
+	   mexErrMsgTxt("unknown boundary detection flags");
+	if(cfg->bc[i])
+	   isbcdet=1;
+     }
+
      if(cfg->medianum){
         for(i=0;i<cfg->medianum;i++)
              if(cfg->prop[i].mus==0.f){
@@ -926,7 +934,7 @@ void mcx_validate_config(Config *cfg){
 		cfg->prop[i].mua*=cfg->unitinmm;
         }
      }
-     if(cfg->issavedet && cfg->detnum==0) 
+     if(cfg->issavedet && cfg->detnum==0 && isbcdet==0)
       	cfg->issavedet=0;
      if(cfg->issavedet==0){
          cfg->issaveexit=0;