bug fix for continuous varying media patch

fangq · May 24, 2019 · 02efc62 · 02efc62
1 parent e1b50b6
commit 02efc62
Show file tree

Hide file tree

Showing 4 changed files with 41 additions and 4 deletions.
diff --git a/mcxlab/mcxlab.m b/mcxlab/mcxlab.m
@@ -351,7 +351,7 @@
 if(nargout>=2)
 
     for i=1:length(varargout{2})
-        if(~isfield(cfg(i),'savedetflag'))
+        if((~isfield(cfg(i),'savedetflag')) || ((isfield(cfg(i),'savedetflag')) && isempty(cfg(i).savedetflag)))
             cfg(i).savedetflag='DP';
             if(isfield(cfg(i),'issaveexit') && cfg(i).issaveexit)
                 cfg(i).savedetflag=[cfg(i).savedetflag,'XV'];

diff --git a/src/mcx_core.cu b/src/mcx_core.cu
@@ -456,9 +456,13 @@ __device__ inline float reflectcoeff(MCXdir *v, float n1, float n2, int flipdir)
 __device__ void updateproperty(Medium *prop, unsigned int mediaid){
 	  if(gcfg->mediaformat<=4)
 	      *((float4*)(prop))=gproperty[mediaid & MED_MASK];
-          else if(gcfg->mediaformat==MEDIA_MUA_FLOAT)
+          else if(gcfg->mediaformat==MEDIA_MUA_FLOAT){
 	      prop->mua=fabs(*((float *)&mediaid));
-	  else if(gcfg->mediaformat==MEDIA_AS_F2H||gcfg->mediaformat==MEDIA_AS_HALF){
+	      if(prop->mua==0.f)
+	          prop->n=1.f;
+	      else
+	          prop->n=gproperty[1].w;
+	  }else if(gcfg->mediaformat==MEDIA_AS_F2H||gcfg->mediaformat==MEDIA_AS_HALF){
 	      union {
                  unsigned int i;
 #if ! defined(__CUDACC_VER_MAJOR__) || __CUDACC_VER_MAJOR__ >= 9
@@ -2235,7 +2239,7 @@ is more than what your have specified (%d), please use the -H option to specify
 	         int j;
 	         for(iter=0;iter<gpu[gpuid].maxgate;iter++)
 		     for(j=0;j<(int)dimlen.z;j++)
-		         mcx_kahanSum(&energyabs[i],&kahanc,cfg->exportfield[iter*dimxyz+(j*cfg->srcnum+i)]*cfg->prop[(uint)cfg->vol[j] & MED_MASK].mua);
+		         mcx_kahanSum(&energyabs[i],&kahanc,cfg->exportfield[iter*dimxyz+(j*cfg->srcnum+i)]*mcx_updatemua((uint)cfg->vol[j],cfg));
 	     }
 	 }
      }

diff --git a/src/mcx_utils.c b/src/mcx_utils.c
@@ -538,6 +538,38 @@ void mcx_normalize(float field[], float scale, int fieldlen, int option, int pid
      *sum=kahant;
  }
 
+ /**
+ * @brief Retrieve mua for different cfg.vol formats to convert fluence back to energy in post-processing 
+ *
+ * @param[out] output: medium absorption coefficient for the current voxel
+ * @param[in] mediaid: medium index of the current voxel
+ * @param[in] cfg: simulation configuration
+ */
+
+ float mcx_updatemua(unsigned int mediaid, Config *cfg){
+     float mua;
+     if(cfg->mediabyte<=4)
+         mua=cfg->prop[mediaid & MED_MASK].mua;
+     else if(cfg->mediabyte==MEDIA_MUA_FLOAT)
+         mua=fabs(*((float *)&mediaid));
+     else if(cfg->mediabyte==MEDIA_ASGN_BYTE){
+         union {
+            unsigned i;
+            unsigned char h[4];
+        } val;
+        val.i=mediaid & MED_MASK;
+        mua=val.h[0]*(1.f/255.f)*(cfg->prop[2].mua-cfg->prop[1].mua)+cfg->prop[1].mua;
+     }else if(cfg->mediabyte==MEDIA_AS_SHORT){
+         union {
+            unsigned int i;
+            unsigned short h[2];
+         } val;
+        val.i=mediaid & MED_MASK;
+        mua=val.h[0]*(1.f/65535.f)*(cfg->prop[2].mua-cfg->prop[1].mua)+cfg->prop[1].mua;
+     }
+     return mua;
+ }
+
 /**
  * @brief Force flush the command line to print the message
  *

diff --git a/src/mcx_utils.h b/src/mcx_utils.h
@@ -257,6 +257,7 @@ int  mcx_isbinstr(const char * str);
 void mcx_progressbar(float percent, Config *cfg);
 void mcx_flush(Config *cfg);
 int  mcx_run_from_json(char *jsonstr);
+float mcx_updatemua(unsigned int mediaid, Config *cfg);
 
 #ifdef MCX_CONTAINER
 #ifdef __cplusplus