add momentum transfer tracking and auto-correlation computation

mmc/trunk/matlab/generate_g1.m

@@ -0,0 +1,89 @@
+function [tau,g1]=generate_g1(fhist,tau, disp_model, DV, lambda, format, varargin)
+%
+%   g1=generate_g1(fhist,tau, disp_model, DV, lambda, format)
+%
+%   Compute simulated electric-field auto-correlation function using
+%   simulated photon pathlengths and scattering momentum transfer
+%
+%   author: Stefan Carp (carp <at> nmr.mgh.harvard.edu)
+%
+%   input:
+%       fhist:      the file name of the output .mch file 
+%       tau:        correlation times at which to compute g1 
+%                   (default: 1e-7 to 1e-1 seconds, log equidistant)
+%       disp_model: displacement model ('brownian', 'random_flow', <custom>)
+%                   (default: brownian, see further explanation below)
+%       disp_var:   value of displacement variable using mm as unit of
+%                   length and s as unit of time
+%                   (default: 1e-7 mm^2/s, see further explanation below)
+%       lambda:     wavelenght of light used in nm
+%                   (default: 785)
+%       format:     the format used to save the .mch file 
+%                   (default: 'float')
+%
+%   output:
+%       g1:         field auto-correlation curves, one for each detector
+%
+%   The displacement model indicates the formula used to compute the root
+%   mean square displacement of scattering particles during a given delay
+%   
+%   brownian:       RMS= 6 * DV * tau; 
+%                   DV(displacement variable)=Db (brownian diffusion coeff)
+%   random_flow:    RMS= DV^2 * tau^2; 
+%                   DV = V (first moment of velocity distribution)
+%   <custom>:       any string other than 'brownian' or 'random_flow' will
+%                   be evaluate as is using Matlab evalf, make sure it uses
+%                   'DV' as the flow related independent variable, tau is
+%                   indexed as tau(J). Any additional parameters can be 
+%                   sent via "varargin"
+%
+%   This file is part of Mesh-Based Monte Carlo
+%   License: GPLv3, see http://mcx.sf.net for details
+%
+
+if nargin<6, format='float'; end
+if nargin<5, lambda=785; end
+if nargin<4, DV=1e-7; end
+if nargin<3, disp_model='brownian'; end
+if nargin<2, tau=logspace(-7,-1,200); end
+
+[mch_data,mch_header]=loadmch(fhist,format);
+temp=strfind(fhist,filesep);
+if isempty(temp), fhist=[pwd filesep fhist]; end
+temp=strfind(fhist,filesep);
+lastslash=temp(end);
+sim_label=fhist(lastslash+1:end-4);
+
+[mua,mus,g,n]=load_mc_prop([fhist(1:lastslash) filesep 'prop_' sim_label '.dat']);
+
+if (mch_header.recordnum-2)~=(2*mch_header.medianum),
+    fprintf('History file does not contain momentum transfer information \n');
+    g1=-1;
+    return;
+end
+
+if strcmp(disp_model,'brownian'),
+    disp_str='rmsdisp=6*DV.*tau(J);';
+elseif strcmp(disp_model,'random_flow'),
+    disp_str='rmsdisp=DV.^2.*tau(J).^2;';
+else 
+    disp_str=['rmsdisp=' disp_model ';'];
+end
+
+k0=2*pi*n/(lambda*1e-6);
+
+g1=zeros(mch_header.detnum,length(tau));
+
+for I=1:mch_header.detnum,
+    idx= find(mch_data(:,1)==I);     
+    fprintf('Processing detector %.0f: %.0f photons\n',I,length(idx));
+
+    for J=1:length(tau),
+        eval(disp_str);
+        g1(I,J)=sum(exp(-(k0.^2.*rmsdisp/3)*mch_data(idx,(3+mch_header.medianum):end)'-mua*mch_data(idx,3:(3+mch_header.medianum-1))'));
+    end
+    g1(I,:)=g1(I,:)./max(g1(I,:));
+end
+
+
+

mmc/trunk/matlab/loadmch.m

@@ -40,14 +40,14 @@
 		end
 		break;
 	end
-	hd=fread(fid,7,'uint');
+	hd=fread(fid,7,'uint'); % version, maxmedia, detnum, colcount, totalphoton, detected, savedphoton
 	if(hd(1)~=1) error('version higher than 1 is not supported'); end
 	unitmm=fread(fid,1,'float32');
 	junk=fread(fid,7,'uint');
 
 	dat=fread(fid,hd(7)*hd(4),format);
 	dat=reshape(dat,[hd(4),hd(7)])';
-	dat(:,3:end)=dat(:,3:end)*unitmm;
+	dat(:,3:(2+hd(2)))=dat(:,3:(2+hd(2)))*unitmm;
 	data=[data;dat];
 	if(isempty(header))
 		header=[hd;unitmm]';

mmc/trunk/src/mcx_utils.c

@@ -32,14 +32,15 @@
 
 const char shortopt[]={'h','E','f','n','t','T','s','a','g','b','B','D',
                  'd','r','S','e','U','R','l','L','I','o','u','C','M',
-		 'i','V','O','\0'};
+                 'i','V','O','m','\0'};
 const char *fullopt[]={"--help","--seed","--input","--photon",
                  "--thread","--blocksize","--session","--array",
                  "--gategroup","--reflect","--reflect3","--debug","--savedet",
                  "--repeat","--save2pt","--minenergy",
                  "--normalize","--skipradius","--log","--listgpu",
                  "--printgpu","--root","--unitinmm","--continuity",
-		 "--method","--interactive","--specular","--outputtype",""};
+                 "--method","--interactive","--specular","--outputtype",
+                 "--momentum",""};
 
 const char debugflag[]={'M','C','B','W','D','I','O','X','A','T','R','P','E','\0'};
 const char raytracing[]={'p','h','b','s','\0'};
@@ -88,7 +89,7 @@ void mcx_initcfg(mcconfig *cfg){
      cfg->srctype=0;
      cfg->isspecular=0;
      cfg->outputtype=otFlux;
-
+     cfg->ismomentum=0;
      cfg->his.version=1;
      cfg->his.unitinmm=1.f;
      memcpy(cfg->his.magic,"MCXH",4);
@@ -465,6 +466,10 @@ void mcx_parsecmd(int argc, char* argv[], mcconfig *cfg){
 		     case 'd':
 		     	        i=mcx_readarg(argc,argv,i,&(cfg->issavedet),"bool");
 		     	        break;
+		     case 'm':
+		                i=mcx_readarg(argc,argv,i,&(cfg->ismomentum),"bool");
+				if (cfg->ismomentum) cfg->issavedet=1;
+				break;			       
 		     case 'C':
 		     	        i=mcx_readarg(argc,argv,i,&(cfg->basisorder),"bool");
 		     	        break;
@@ -566,6 +571,7 @@ where possible parameters include (the first item in [] is the default value)\n\
  -e [0.|float] (--minenergy)   minimum energy level to trigger Russian roulette\n\
  -U [1|0]      (--normalize)   1 to normalize the fluence to unitary,0 save raw\n\
  -d [0|1]      (--savedet)     1 to save photon info at detectors,0 not to save\n\
+ -m [0|1]      (--momentum)    1 to save photon momentum transfer,0 not to save\n\
  -S [1|0]      (--save2pt)     1 to save the fluence field, 0 do not save\n\
  -C [1|0]      (--basisorder)  1 piece-wise-linear basis for fluence,0 constant\n\
  -V [0|1]      (--specular)    1 source located in the background,0 inside mesh\n\

mmc/trunk/src/mcx_utils.h

@@ -132,6 +132,7 @@ typedef struct MMC_config{
         char isref3;        /**<1 considering maximum 3 ref. interfaces; 0 max 2 ref*/
 	char isnormalized;  /**<1 to normalize the fluence, 0 for raw fluence*/
 	char issavedet;     /**<1 to count all photons hits the detectors*/
+	char ismomentum;    /**<1 to save momentum transfer for detected photons, implies issavedet=1*/
 	char issave2pt;     /**<1 to save the 2-point distribution, 0 do not save*/
 	char isgpuinfo;     /**<1 to print gpu info when attach, 0 do not print*/
 	char isspecular;    /**<1 calculate the initial specular ref if outside the mesh, 0 do not calculate*/

mmc/trunk/src/simpmesh.c

@@ -428,7 +428,9 @@ void tracer_clear(raytracer *tracer){
 	tracer->mesh=NULL;
 }
 
-float mc_next_scatter(float g, float3 *dir,RandType *ran, RandType *ran0, mcconfig *cfg){
+
+float mc_next_scatter(float g, float3 *dir,RandType *ran, RandType *ran0, mcconfig *cfg, float *pmom){
+
     float nextslen;
     float sphi,cphi,tmp0,theta,stheta,ctheta,tmp1;
     float3 p;
@@ -484,6 +486,9 @@ float mc_next_scatter(float g, float3 *dir,RandType *ran, RandType *ran0, mcconf
 	p.y=stheta*sphi;
 	p.z=(dir->z>0.f)?ctheta:-ctheta;
     }
+    if (cfg->ismomentum)
+        pmom[0]+=(1.f-ctheta);
+
     dir->x=p.x;
     dir->y=p.y;
     dir->z=p.z;
@@ -536,7 +541,7 @@ void mesh_savedetphoton(float *ppath, int count, mcconfig *cfg){
         cfg->his.detected=count;
 	cfg->his.savedphoton=count;
 	cfg->his.detnum=cfg->detnum;
-        cfg->his.colcount=cfg->his.maxmedia+2; /*column count=maxmedia+2*/
+        cfg->his.colcount=(1+(cfg->ismomentum>0))*cfg->his.maxmedia+2; /*column count=maxmedia+2*/
 
 	fwrite(&(cfg->his),sizeof(history),1,fp);
 	fwrite(ppath,sizeof(float),count*cfg->his.colcount,fp);

mmc/trunk/src/simpmesh.h

@@ -113,8 +113,8 @@ void tracer_init(raytracer *tracer,tetmesh *mesh,char methodid);
 void tracer_build(raytracer *tracer);
 void tracer_prep(raytracer *tracer,mcconfig *cfg);
 void tracer_clear(raytracer *tracer);
-float mc_next_scatter(float g, float3 *dir,RandType *ran,RandType *ran0,mcconfig *cfg);
 
+float mc_next_scatter(float g, float3 *dir,RandType *ran,RandType *ran0,mcconfig *cfg,float *pmom);
 
 static inline void vec_add(float3 *a,float3 *b,float3 *res){
 	res->x=a->x+b->x;

mmc/trunk/src/tetray.c

@@ -85,8 +85,9 @@ int main(int argc, char**argv){
 #pragma omp parallel private(ran0,ran1,threadid)
 {
 	visitor visit={0.f,1.f/cfg.tstep,DET_PHOTON_BUF,0,NULL};
+	int buflen=(1+((cfg.ismomentum)>0))*mesh.prop+2;
 	if(cfg.issavedet) 
-	    visit.partialpath=(float*)calloc(visit.detcount*(mesh.prop+2),sizeof(float));
+	    visit.partialpath=(float*)calloc(visit.detcount*buflen,sizeof(float));
 #ifdef _OPENMP
 	threadid=omp_get_thread_num();	
 #endif
@@ -109,12 +110,12 @@ int main(int argc, char**argv){
 		master.detcount+=visit.bufpos;
             #pragma omp barrier
 	    if(threadid==0)
-	        master.partialpath=(float*)calloc(master.detcount*(mesh.prop+2),sizeof(float));
+		master.partialpath=(float*)calloc(master.detcount*buflen,sizeof(float));
             #pragma omp barrier
             #pragma omp critical
             {
-		memcpy(master.partialpath+master.bufpos*(mesh.prop+2),
-		    visit.partialpath,visit.bufpos*(mesh.prop+2)*sizeof(float));
+		memcpy(master.partialpath+master.bufpos*buflen,
+		       visit.partialpath,visit.bufpos*buflen*sizeof(float));
 		master.bufpos+=visit.bufpos;
             }
             #pragma omp barrier

mmc/trunk/src/tettracing.c

@@ -661,7 +661,7 @@ float onephoton(int id,raytracer *tracer,tetmesh *mesh,mcconfig *cfg,
 			r.vec.z=cosf(aangle);
 		}
 	}
-	r.partialpath=(float*)calloc(mesh->prop+1,sizeof(float));
+	r.partialpath=(float*)calloc((1+(cfg->ismomentum>0))*mesh->prop+1,sizeof(float));
 
 	tracercore=engines[0];
 	if(cfg->method>=0 && cfg->method<4)
@@ -763,18 +763,20 @@ float onephoton(int id,raytracer *tracer,tetmesh *mesh,mcconfig *cfg,
 		else
 			break;
 	    }
-	    r.slen=mc_next_scatter(mesh->med[mesh->type[r.eid-1]].g,&r.vec,ran,ran0,cfg);
+	    r.slen=mc_next_scatter(mesh->med[mesh->type[r.eid-1]].g,&r.vec,ran,ran0,cfg,
+	                           r.partialpath+mesh->prop+mesh->type[r.eid-1]);
             r.partialpath[0]++;
 	}
 	if(cfg->issavedet && exitdet>0){
-		int offset=visit->bufpos*(mesh->prop+2);
+		int buflen=(1+(cfg->ismomentum>0))*mesh->prop+2;
+		int offset=visit->bufpos*buflen;
 		if(visit->bufpos>=visit->detcount){
 		    visit->detcount+=DET_PHOTON_BUF;
 		    visit->partialpath=(float *)realloc(visit->partialpath,
-		               visit->detcount*(mesh->prop+2)*sizeof(float));
+				visit->detcount*buflen*sizeof(float));
 		}
 		visit->partialpath[offset]=exitdet;
-		memcpy(visit->partialpath+offset+1,r.partialpath, (mesh->prop+1)*sizeof(float));
+	        memcpy(visit->partialpath+offset+1,r.partialpath,(buflen-1)*sizeof(float));
 		visit->bufpos++;
 	}
 	free(r.partialpath);