initial support for hybrid optical properties: tissue type label comb…

…ined with continous optical properties

mcxlab/examples/demo_label_continous_hybrid.m

@@ -0,0 +1,153 @@
+clear
+clc
+
+%% case 1: homogenous media, update the refractive index from n_old to n_new for all voxels
+clear cfg1
+n_old=1.37; 
+    n_new=1.55;
+
+cfg1.nphoton=1e8;
+cfg1.srcpos=[30 30 1];
+cfg1.srcdir=[0 0 1];
+cfg1.gpuid=1;
+cfg1.autopilot=1;
+cfg1.tstart=0;
+cfg1.tend=5e-9;
+cfg1.tstep=5e-9;
+cfg1.isreflect=1;
+
+% conventional way to define optical properties: label each voxel with
+% tissue types
+cfg1.vol=uint8(ones(60,60,60));
+    cfg1.prop=[0 0 1 1;0.005 1 0 n_new];
+flux1_conventional=mcxlab(cfg1);
+
+% tissue-type based continuous varying media
+cfg1.prop=[0 0 1 1;0.005 1 0 n_old];                     % n=1.37 will be updated during simulation
+    property_substitute=single(4*ones(size(cfg1.vol)));  % 0-no update, 1-update mua, 2-update mus, 3-update g, 4-update n
+    label=single(ones(size(cfg1.vol)));                  % label-based tissue type 0-255
+    replace_property=single(n_new*ones(size(cfg1.vol))); % values to substitute corresponding optical properties
+    cfg1.vol=permute(cat(4,replace_property,label,property_substitute),[4 1 2 3]); % Concatenated as a 4-D single array 
+flux1_new=mcxlab(cfg1);
+
+% compare fluence maps between two approaches
+figure(1);
+clines=-7:1:-1;
+contourf(log10(squeeze(flux1_conventional.data(31,:,:))*cfg1.tstep),clines,'r-');axis equal;
+hold on;
+contour(log10(squeeze(flux1_new.data(31,:,:))*cfg1.tstep),clines,'w--');
+legend('conventional','new');
+colorbar;
+title("homogeneous cube");
+
+%% case 2: Two-layer slab: The refractive indices of tissue type 1 and 2 will be updated to 1.45 and 1.6 respectively
+clear cfg2
+n_old_1=1.37; 
+    n_new_1=1.45;
+
+n_old=1.37; 
+    n_new_2=1.6;
+
+cfg2.nphoton=1e8;
+cfg2.srcpos=[30 30 1];
+cfg2.srcdir=[0 0 1];
+cfg2.gpuid=1;
+cfg2.autopilot=1;
+cfg2.tstart=0;
+cfg2.tend=5e-9;
+cfg2.tstep=5e-9;
+cfg2.isreflect=1;
+
+% conventional way to define optical properties: label each voxel with
+% tissue types
+cfg2.vol=uint8(ones(60,60,60));
+    cfg2.vol(:,:,31:60)=2;
+cfg2.prop=[0 0 1 1;0.005 1 0 n_new_1;0.01 2 0 n_new_2];
+flux2_conventional=mcxlab(cfg2);
+
+% tissue-type based continuous varying media
+cfg2.prop=[0 0 1 1;0.005 1 0 n_old_1;0.01 2 0 n_old];      % n=1.37 will be updated during the simulation
+    property_substitute=single(4*ones(size(cfg2.vol)));    % 0-no update, 1-update mua, 2-update mus, 3-update g, 4-update n
+    label=single(ones(size(cfg2.vol)));                    % label-based tissue type 0-255
+        label(:,:,31:60)=2;
+    replace_property=single(n_new_1*ones(size(cfg2.vol))); % values to substitute corresponding optical properties
+        replace_property(:,:,31:60)=n_new_2;
+    cfg2.vol=permute(cat(4,replace_property,label,property_substitute),[4 1 2 3]); % Concatenated as a 4-D array 
+flux2_new=mcxlab(cfg2);
+
+clines=[-10:0.5:-2];
+% compare
+figure(2);
+contourf(log10(squeeze(flux2_conventional.data(31,:,:))*cfg2.tstep),clines,'r-');axis equal;
+hold on;
+contour(log10(squeeze(flux2_new.data(31,:,:))*cfg2.tstep),clines,'w--');
+legend('conventional','new');
+title("two layer slab");
+
+%% case 3: Two-layer slab: mua of tissue type 1 will be updated to 0.05 while the refractive indices of tissue type 2 will be updated to [1.30,1.35,1.40,1.45,1.50,1.55]
+clear cfg3
+mua_old_1=0.005; 
+    mua_new_1=0.05;
+
+n_old=1.37;
+    n_new_1=1.30;
+    n_new_2=1.35;
+    n_new_3=1.40;
+    n_new_4=1.45;
+    n_new_5=1.50;
+    n_new_6=1.55;
+
+cfg3.nphoton=1e8;
+cfg3.srcpos=[30 30 1];
+cfg3.srcdir=[0 0 1];
+cfg3.gpuid=1;
+cfg3.autopilot=1;
+cfg3.tstart=0;
+cfg3.tend=5e-9;
+cfg3.tstep=5e-9;
+cfg3.isreflect=1;
+
+% conventional way to define optical properties: label each voxel with
+% tissue types
+cfg3.vol=uint8(ones(60,60,60));
+    for i=1:6     % each layer
+        cfg3.vol(:,:,30+i*5)=i+1;
+        cfg3.vol(:,:,30+i*5-1)=i+1;
+        cfg3.vol(:,:,30+i*5-2)=i+1;
+        cfg3.vol(:,:,30+i*5-3)=i+1;
+        cfg3.vol(:,:,30+i*5-4)=i+1;
+    end
+cfg3.prop=[0 0 1 1;
+           mua_new_1 1 0 1.37;
+           0.01 2 0 n_new_1;
+           0.01 2 0 n_new_2;
+           0.01 2 0 n_new_3;
+           0.01 2 0 n_new_4;
+           0.01 2 0 n_new_5;
+           0.01 2 0 n_new_6];
+flux2_conventional=mcxlab(cfg3);
+
+% tissue-type based continuous varying media
+cfg3.prop=[0 0 1 1;0.005 1 0 1.37;0.01 2 0 n_old];         % n=1.37 will be updated during the simulation
+    property_substitute=single(ones(size(cfg3.vol)));      % 0-no update, 1-update mua, 2-update mus, 3-update g, 4-update n
+        property_substitute(:,:,31:60)=4;
+    label=single(ones(size(cfg3.vol)));                    % label-based tissue type 0-255
+        label(:,:,31:60)=2;
+    replace_property=single(mua_new_1*ones(size(cfg3.vol))); % values to substitute corresponding optical properties
+        replace_property(:,:,31:35)=n_new_1;
+        replace_property(:,:,36:40)=n_new_2;
+        replace_property(:,:,41:45)=n_new_3;
+        replace_property(:,:,46:50)=n_new_4;
+        replace_property(:,:,51:55)=n_new_5;
+        replace_property(:,:,56:60)=n_new_6;
+    cfg3.vol=permute(cat(4,replace_property,label,property_substitute),[4 1 2 3]); % Concatenated as a 4-D array 
+flux2_new=mcxlab(cfg3);
+
+clines=[-12:1:-2];
+% compare
+figure(3);
+contourf(log10(squeeze(flux2_conventional.data(31,:,:))*cfg3.tstep),clines,'r-');axis equal;
+hold on;
+contour(log10(squeeze(flux2_new.data(31,:,:))*cfg3.tstep),clines,'w--');
+legend('conventional','new');
+title("multi-layered slab");

src/mcx_const.h

@@ -49,6 +49,7 @@
 #define MCX_DEBUG_MOVE      2
 #define MCX_DEBUG_PROGRESS  4
 
+#define MEDIA_LABEL_HALF      99   /**<  media format: {[float32: 1/2/3/4][float32: type][float32: mua/mus/g/n]} -> 32bit:{[byte: 1/2/3/4][byte: type][half: mua/mus/g/n]} */
 #define MEDIA_AS_F2H          100  /**<  media format: {[float32: mua][float32: mus]} -> 32bit:{[half: mua],{half: mus}} */
 #define MEDIA_MUA_FLOAT       101  /**<  media format: 32bit:{[float32: mua]} */
 #define MEDIA_AS_HALF         102  /**<  media format: 32bit:{[half: mua],[half: mus]} */

src/mcx_core.cu

@@ -463,9 +463,25 @@ __device__ inline float reflectcoeff(MCXdir *v, float n1, float n2, int flipdir)
  */
 
 __device__ void updateproperty(Medium *prop, unsigned int mediaid){
-	  if(gcfg->mediaformat<=4)
+	  if(gcfg->mediaformat<=4){
 	      *((float4*)(prop))=gproperty[mediaid & MED_MASK];
-          else if(gcfg->mediaformat==MEDIA_MUA_FLOAT){
+	  }else if(gcfg->mediaformat==MEDIA_LABEL_HALF){
+	      union{
+	         unsigned int i;
+#if ! defined(__CUDACC_VER_MAJOR__) || __CUDACC_VER_MAJOR__ >= 9
+                 __half_raw h[2];
+#else
+                 half h[2];
+#endif
+		 unsigned char c[4];	  
+	      } val;
+	      val.i=mediaid & MED_MASK;
+	      *((float4*)(prop))=gproperty[(unsigned int)(val.c[2])];
+	      if((unsigned int)(val.c[3])>0 && (unsigned int)(val.c[3])<5){
+	         float *p=(float*)(prop);
+	         p[(unsigned int)(val.c[3])-1]=fabs(__half2float(val.h[0]));
+	      }
+          }else if(gcfg->mediaformat==MEDIA_MUA_FLOAT){
 	      prop->mua=fabs(*((float *)&mediaid));
               prop->n=gproperty[!(mediaid & MED_MASK)==0].w;
 	  }else if(gcfg->mediaformat==MEDIA_AS_F2H||gcfg->mediaformat==MEDIA_AS_HALF){
@@ -1431,7 +1447,9 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
 
           /** do boundary reflection/transmission */
 	  if(isreflect){
-	      if(((gcfg->doreflect && (isdet & 0xF)==0) || (isdet & 0x1)) && n1!=gproperty[(mediaid>0 && gcfg->mediaformat>=100)?1:mediaid].w){
+	      if(gcfg->mediaformat==MEDIA_LABEL_HALF)
+	          updateproperty(&prop,mediaid); ///< optical property across the interface
+	      if(((gcfg->doreflect && (isdet & 0xF)==0) || (isdet & 0x1)) && n1!=((gcfg->mediaformat==MEDIA_LABEL_HALF)? (prop.n):(gproperty[(mediaid>0 && gcfg->mediaformat>=100)?1:mediaid].w))){
 	          float Rtotal=1.f;
 	          float cphi,sphi,stheta,ctheta,tmp0,tmp1;
 
@@ -1453,7 +1471,7 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
        	       		Rtotal=(Rtotal+(ctheta-stheta)/(ctheta+stheta))*0.5f;
 	        	GPUDEBUG(("Rtotal=%f\n",Rtotal));
                   } ///< else, total internal reflection
-	          if(Rtotal<1.f && (((isdet & 0xF)==0 && gproperty[mediaid].w>=1.f) || isdet==bcReflect) && rand_next_reflect(t)>Rtotal){ // do transmission
+	          if(Rtotal<1.f && (((isdet & 0xF)==0 && ((gcfg->mediaformat==MEDIA_LABEL_HALF) ? prop.n:gproperty[mediaid].w) >= 1.f) || isdet==bcReflect) && rand_next_reflect(t)>Rtotal){ // do transmission
                         transmit(&v,n1,prop.n,flipdir);
                         if(mediaid==0){ // transmission to external boundary
                             GPUDEBUG(("transmit to air, relaunch\n"));
@@ -1482,7 +1500,8 @@ kernel void mcx_main_loop(uint media[],OutputType field[],float genergy[],uint n
         	  	updateproperty(&prop,mediaid); ///< optical property across the interface
                   	n1=prop.n;
 		  }
-	      }
+	      }else if(gcfg->mediaformat==MEDIA_LABEL_HALF)
+	          updateproperty(&prop,mediaidold); ///< optical property across the interface
 	  }
      }
 

src/mcxlab.cpp

@@ -480,6 +480,8 @@ void mcx_set_field(const mxArray *root,const mxArray *item,int idx, Config *cfg)
 	}else if(mxGetNumberOfDimensions(item)==4){ // if dimension is 4D, 1st dim is the property records: mua/mus/g/n
 	    if((mxIsUint8(item) || mxIsInt8(item)) && arraydim[0]==4) // if 4D byte array has a 1st dim of 4
 		 cfg->mediabyte=MEDIA_ASGN_BYTE;
+	    else if(arraydim[0]==3)
+		 cfg->mediabyte=MEDIA_LABEL_HALF;
 	    else if((mxIsUint16(item) || mxIsInt16(item)) && arraydim[0]==2)// if 4D short array has a 1st dim of 2
 		 cfg->mediabyte=MEDIA_AS_SHORT;
 	    else if(mxIsSingle(item) && arraydim[0]==2) // if 4D float32 array has a 1st dim of 2
@@ -541,6 +543,31 @@ void mcx_set_field(const mxArray *root,const mxArray *item,int idx, Config *cfg)
 
 	            cfg->vol[i]=f2h.i[0];
 		}
+	    }else if(cfg->mediabyte==MEDIA_LABEL_HALF){
+		float *val=(float *)mxGetPr(item);
+		union{
+		    float f[3];
+		    unsigned int i[3];
+		    unsigned short h[2];
+		    unsigned char c[4];
+		} f2bh;
+		unsigned short tmp;
+		for(i=0;i<dimxyz;i++){
+		    f2bh.f[0]=val[i*3];
+		    f2bh.f[1]=val[i*3+1];
+		    f2bh.f[2]=val[i*3+2];
+
+		    f2bh.h[0] = (f2bh.i[0] >> 31) << 5;
+		    tmp = (f2bh.i[0] >> 23) & 0xff;
+		    tmp = (tmp - 0x70) & ((unsigned int)((int)(0x70 - tmp) >> 4) >> 27);
+		    f2bh.h[0] = (f2bh.h[0] | tmp) << 10;
+		    f2bh.h[0] |= (f2bh.i[0] >> 13) & 0x3ff;
+
+		    f2bh.c[3]=(unsigned char)(f2bh.f[2]);
+		    f2bh.c[2]=(unsigned char)(f2bh.f[1]);
+
+	            cfg->vol[i]=f2bh.i[0];
+		}
 	    }
 	}
         printf("mcx.dim=[%d %d %d];\n",cfg->dim.x,cfg->dim.y,cfg->dim.z);