feat: output prop. along with det. photon profile

After polarized MCXLAB simulations, return mus and g of each tissue
type to the 'prop' field of detphoton for post-processing.

mcxlab/mcxlab.m

@@ -466,12 +466,9 @@
                 flags{end+1}=cfg(i).srcnum;
             end
             newdetp=mcxdetphoton(detp,medianum,flags{:});
-            if(isfield(cfg(i),'polprop') && ~isempty(cfg(i).polprop))
-                newdetp.prop=zeros(medianum+1,4);
-                newdetp.prop(1,:)=cfg(i).prop(1,:);
-                newdetp.prop(2:end,1)=cfg(i).polprop(:,1);
-            else
-                newdetp.prop=cfg(i).prop;
+            newdetp.prop=cfg(i).prop;
+            if(isfield(cfg(i),'polprop') && ~isempty(cfg(i).polprop)) && isfield(varargout{1}(i),'prop')
+                newdetp.prop(2:end,:)=varargout{1}(i).prop(:,2:end)';
             end
             if(isfield(cfg(i),'unitinmm'))
                 newdetp.unitinmm=cfg(i).unitinmm;

src/mcx_utils.c

@@ -1373,8 +1373,6 @@ void mcx_prep_polarized(Config *cfg){
     cfg->smatrix=(float4 *)malloc(cfg->polmedianum*NANGLES*sizeof(float4));
     Medium *prop=cfg->prop;
     POLMedium *polprop=cfg->polprop;
-    FILE *pfile;
-    pfile=fopen("Mie_outputs.dat","w");
 
     for(int i=0;i<cfg->polmedianum;i++){
         prop[i+1].mua=polprop[i].mua;
@@ -1402,16 +1400,10 @@ void mcx_prep_polarized(Config *cfg){
         /* g will store the index that points to the corresponding smatrix */
         prop[i+1].g=(float)i;
 
-        /* save Mie function outputs to a file */
-        if(pfile!=NULL){
-            fprintf(pfile,"Tissue type %d: ",i+1);
-            fprintf(pfile,"mua=%5.5f mm^(-1), radius=%5.5f micron, rho=%5.5f micron^(-3), n_sph=%5.5f, n_med=%5.5f\n",
-                polprop[i].mua,polprop[i].r,polprop[i].rho,polprop[i].nsph,polprop[i].nmed);
-            fprintf(pfile,"Mie function outputs: mus=%5.5f mm^(-1), g=%5.5f\n\n",prop[i+1].mus,g);
-        }
+        /* store g in polprop, which will be output to detphoton later for post-processing */
+        polprop[i].mua=g;
     }
     free(mu);
-    if(pfile) fclose(pfile);
 }
 
 /**

src/mcxlab.cpp

@@ -102,7 +102,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
   int        errorflag=0;
   int        threadid=0;
   const char       *outputtag[]={"data"};
-  const char       *datastruct[]={"data","stat","dref"};
+  const char       *datastruct[]={"data","stat","dref","prop"};
   const char       *statstruct[]={"runtime","nphoton","energytot","energyabs","normalizer","unitinmm","workload"};
   const char       *gpuinfotag[]={"name","id","devcount","major","minor","globalmem",
                                   "constmem","sharedmem","regcount","clock","sm","core",
@@ -171,7 +171,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
    * The function can return 1-5 outputs (i.e. the LHS)
    */
   if(nlhs>=1)
-      plhs[0] = mxCreateStructMatrix(ncfg,1,3,datastruct);
+      plhs[0] = mxCreateStructMatrix(ncfg,1,4,datastruct);
   if(nlhs>=2)
       plhs[1] = mxCreateStructMatrix(ncfg,1,1,outputtag);
   if(nlhs>=3)
@@ -386,6 +386,18 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
             mxSetFieldByNumber(stat,0,6, val);
 
 	    mxSetFieldByNumber(plhs[0],jstruct,1, stat);
+
+            /** return the final optical properties for polarized MCX simulation */
+            if(cfg.polprop){
+                for(int i=0;i<cfg.polmedianum;i++){
+                    // remember that g was replaced by the index of Mueller matrix, now restore it
+                    cfg.prop[i+1].g=cfg.polprop[i].mua;
+                }
+
+                dimtype propdim[2]={4,cfg.medianum};
+                mxSetFieldByNumber(plhs[0],jstruct,3, mxCreateNumericArray(2,propdim,mxSINGLE_CLASS,mxREAL));
+                memcpy((float*)mxGetPr(mxGetFieldByNumber(plhs[0],jstruct,3)),cfg.prop,cfg.medianum*4*sizeof(float));
+            }
         }
     }catch(const char *err){
       mexPrintf("Error: %s\n",err);