prepare for axis-based bounary condition

src/mcx_core.cu

@@ -1745,7 +1745,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),sizeof(uint2));
      Vvox=cfg->steps.x*cfg->steps.y*cfg->steps.z; /*Vvox: voxel volume in mm^3*/
 
      if(cfg->seed>0)

src/mcx_core.h

@@ -147,6 +147,7 @@ typedef struct  __align__(16) KernelParams {
   unsigned int gscatter;             /**< how many scattering events after which mus/g can be approximated by mus' */
   unsigned int is2d;                 /**< is the domain a 2D slice? */
   int replaydet;                     /**< select which detector to replay, 0 for all, -1 save all separately */
+  unsigned char bc[2][3];              /**< is the domain a 2D slice? */
 }MCXParam;
 
 void mcx_run_simulation(Config *cfg,GPUInfo *gpu);

src/mcx_utils.c

@@ -74,9 +74,9 @@
  * Array terminates with '\0'.
  */
 
-const char shortopt[]={'h','i','f','n','t','T','s','a','g','b','B','z','u','H','P','N',
+const char shortopt[]={'h','i','f','n','t','T','s','a','g','b','-','z','u','H','P','N',
                  'd','r','S','p','e','U','R','l','L','-','I','-','G','M','A','E','v','D',
-		 'k','q','Y','O','F','-','-','x','X','-','-','m','V','\0'};
+		 'k','q','Y','O','F','-','-','x','X','-','-','m','V','B','\0'};
 
 /**
  * Long command line options
@@ -93,7 +93,7 @@ const char *fullopt[]={"--help","--interactive","--input","--photon",
 		 "--seed","--version","--debug","--voidtime","--saveseed",
 		 "--replaydet","--outputtype","--outputformat","--maxjumpdebug",
                  "--maxvoidstep","--saveexit","--saveref","--gscatter","--mediabyte",
-                 "--momentum","--specular",""};
+                 "--momentum","--specular","--bc",""};
 
 /**
  * Output data types
@@ -125,6 +125,16 @@ const char debugflag[]={'R','M','P','\0'};
 
 const char *outputformat[]={"mc2","nii","hdr","ubj",""};
 
+/**
+ * Boundary condition (BC) types
+ * r: Fresnel boundary
+ * a: total absorption BC
+ * m: total reflection (mirror) BC
+ * c: cylic BC
+ */
+
+const char boundarycond[]={'r','a','m','c','\0'};
+
 /**
  * Source type specifier
  * User can specify the source type using a string
@@ -221,6 +231,7 @@ void mcx_initcfg(Config *cfg){
      cfg->isspecular=0;
      cfg->dx=cfg->dy=cfg->dz=NULL;
      cfg->gscatter=1e9;     /** by default, honor anisotropy for all scattering, use --gscatter to reduce it */
+     memset(&(cfg->bc),0,sizeof(uint2));
      memset(&(cfg->srcparam1),0,sizeof(float4));
      memset(&(cfg->srcparam2),0,sizeof(float4));
      memset(cfg->deviceid,0,MAX_DEVICE);
@@ -644,6 +655,7 @@ void mcx_writeconfig(char *fname, Config *cfg){
  */
 
 void mcx_prepdomain(char *filename, Config *cfg){
+     char *bc;
      if(filename[0] || cfg->vol){
         if(cfg->vol==NULL){
 	     mcx_loadvolume(filename,cfg);
@@ -694,6 +706,11 @@ void mcx_prepdomain(char *filename, Config *cfg){
      for(int i=0;i<MAX_DEVICE;i++)
         if(cfg->deviceid[i]=='0')
            cfg->deviceid[i]='\0';
+
+     bc=(char*)(&cfg->bc);
+     for(int i=0;i<6;i++)
+        if(bc[i]>='A' && mcx_lookupindex(bc+i,boundarycond))
+	   MCX_ERROR(-4,"unknown boundary condition specifier");
 }
 
 /**
@@ -1755,9 +1772,13 @@ void mcx_parsecmd(int argc, char* argv[], Config *cfg){
 		     case 'b':
 		     	        i=mcx_readarg(argc,argv,i,&(cfg->isreflect),"char");
 				cfg->isref3=cfg->isreflect;
+				if(cfg->isreflect)
+				    memset(&(cfg->bc),bcReflect,sizeof(uint2));
+				else
+				    memset(&(cfg->bc),bcAbsorb,sizeof(uint2));
 		     	        break;
                      case 'B':
-                                i=mcx_readarg(argc,argv,i,&(cfg->isrefint),"char");
+                                i=mcx_readarg(argc,argv,i,&(cfg->bc),"string");
                                	break;
 		     case 'd':
 		     	        i=mcx_readarg(argc,argv,i,&(cfg->issavedet),"char");
@@ -1892,6 +1913,8 @@ void mcx_parsecmd(int argc, char* argv[], Config *cfg){
                                      i=mcx_readarg(argc,argv,i,&(cfg->faststep),"char");
                                 else if(strcmp(argv[i]+2,"root")==0)
                                      i=mcx_readarg(argc,argv,i,cfg->rootpath,"string");
+                                else if(strcmp(argv[i]+2,"reflectin")==0)
+                                     i=mcx_readarg(argc,argv,i,&(cfg->isrefint),"char");
                                 else
                                      MCX_FPRINTF(cfg->flog,"unknown verbose option: --%s\n",argv[i]+2);
 		     	        break;

src/mcx_utils.h

@@ -51,6 +51,7 @@
 enum TOutputType {otFlux, otFluence, otEnergy, otJacobian, otWP, otDCS};   /**< types of output */
 enum TMCXParent  {mpStandalone, mpMATLAB};                          /**< whether MCX is run in binary or mex mode */
 enum TOutputFormat {ofMC2, ofNifti, ofAnalyze, ofUBJSON};           /**< output data format */
+enum TBoundary {bcReflect, bcAbsorb, bcMirror, bcCylic};            /**< boundary conditions */
 
 /**
  * The structure to store optical properties
@@ -210,6 +211,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 */
+	uint2 bc;                    /**<boundary condition flag*/
 } Config;
 
 #ifdef	__cplusplus

src/mcxlab.cpp

@@ -601,6 +601,13 @@ void mcx_set_field(const mxArray *root,const mxArray *item,int idx, Config *cfg)
         jsonshapes=new char[len+1];
         mxGetString(item, jsonshapes, len+1);
         jsonshapes[len]='\0';
+    }else if(strcmp(name,"bc")==0){
+        int len=mxGetNumberOfElements(item);
+        if(!mxIsChar(item) || len==0 || len>7)
+             mexErrMsgTxt("the 'bc' field must be a non-empty string");
+
+        mxGetString(item, (char *)&(cfg->bc), len+1);
+        ((char *)&(cfg->bc))[len]='\0';
     }else if(strcmp(name,"detphotons")==0){
         arraydim=mxGetDimensions(item);
 	dimdetps[0]=arraydim[0];
@@ -731,6 +738,7 @@ void mcx_replay_prep(Config *cfg){
 
 void mcx_validate_config(Config *cfg){
      int i,gates,idx1d;
+     const char boundarycond[]={'r','a','m','c','\0'};
 
      if(!cfg->issrcfrom0){
         cfg->srcpos.x--;cfg->srcpos.y--;cfg->srcpos.z--; /*convert to C index, grid center*/
@@ -780,8 +788,13 @@ void mcx_validate_config(Config *cfg){
      if(cfg->replaydet==-1 && cfg->detnum==1)
         cfg->replaydet=1;
 
+     char *bc=(char*)(&cfg->bc);
+     for(i=0;i<6;i++)
+        if(bc[i]>='A' && mcx_lookupindex(bc+i,boundarycond))
+	   mexErrMsgTxt("unknown boundary condition specifier");
+
      if(cfg->medianum){
-        for(int i=0;i<cfg->medianum;i++)
+        for(i=0;i<cfg->medianum;i++)
              if(cfg->prop[i].mus==0.f){
 	         cfg->prop[i].mus=EPS;
 		 cfg->prop[i].g=1.f;