diff --git a/README.md b/README.md index 7b8c6103..40359136 100644 --- a/README.md +++ b/README.md @@ -1337,7 +1337,7 @@ to adpot the newly implemented JSON/.jdat format for easy data sharing. The .mch file contains a 256 byte binary header, followed by a 2-D numerical array of dimensions `#savedphoton * #colcount` as recorded in the header. ``` - typedef struct MCXHistoryHeader{ +typedef struct MCXHistoryHeader{ char magic[4]; // magic bits= 'M','C','X','H' unsigned int version; // version of the mch file format unsigned int maxmedia; // number of media in the simulation @@ -1346,14 +1346,15 @@ of dimensions `#savedphoton * #colcount` as recorded in the header. unsigned int totalphoton; // how many total photon simulated unsigned int detected; // how many photons are detected (not necessarily all saved) unsigned int savedphoton; // how many detected photons are saved in this file - float unitinmm; // what is the voxel size of the simulation + float unitinmm; // what is the voxel size of the simulation unsigned int seedbyte; // how many bytes per RNG seed - float normalizer; // what is the normalization factor + float normalizer; // what is the normalization factor int respin; // if positive, repeat count so total photon=totalphoton*respin; if negative, total number is processed in respin subset unsigned int srcnum; // number of sources for simultaneous pattern sources unsigned int savedetflag; // number of sources for simultaneous pattern sources - int reserved[2]; // reserved fields for future extension - } History; + unsigned int totalsource; // total source number when multiple sources are defined + int reserved[1]; // reserved fields for future extension +} History; ``` When the `-q` flag is set to 1, the detected photon initial seeds are also stored following the detected photon data, consisting of a 2-D byte array of `#savedphoton * #seedbyte`. @@ -1369,36 +1370,36 @@ When `-F jnii` is specified, instead of saving the detected photon into the lega a .jdat file is written, which is a pure JSON file. This file contains a hierachical data record of the following JSON structure ```` - { - "MCXData": { +{ + "MCXData":{ "Info":{ "Version": - "MediaNum": - "DetNum": - ... - "Media":{ - ... - } + "MediaNum": + "DetNum": + ... + "Media":{ + ... + } }, "PhotonData":{ "detid": - "nscat": - "ppath": - "mom": - "p": - "v": - "w0": + "nscat": + "ppath": + "mom": + "p": + "v": + "w0": }, "Trajectory":{ "photonid": - "p": - "w0": + "p": + "w0": }, "Seed":[ ... ] } - } +} ```` where "Info" is required, and other subfields are optional depends on users' input. Each subfield in this file may contain JData 1-D or 2-D array constructs to allow diff --git a/README.txt b/README.txt index 7e48af46..ae8c9aad 100644 --- a/README.txt +++ b/README.txt @@ -1328,21 +1328,22 @@ The .mch file contains a 256 byte binary header, followed by a 2-D numerical arr of dimensions #savedphoton * #colcount as recorded in the header. typedef struct MCXHistoryHeader{ - char magic[4]; /**< magic bits= 'M','C','X','H' */ - unsigned int version; /**< version of the mch file format */ - unsigned int maxmedia; /**< number of media in the simulation */ - unsigned int detnum; /**< number of detectors in the simulation */ - unsigned int colcount; /**< how many output files per detected photon */ - unsigned int totalphoton; /**< how many total photon simulated */ - unsigned int detected; /**< how many photons are detected (not necessarily all saved) */ - unsigned int savedphoton; /**< how many detected photons are saved in this file */ - float unitinmm; /**< what is the voxel size of the simulation */ - unsigned int seedbyte; /**< how many bytes per RNG seed */ - float normalizer; /**< what is the normalization factor */ - int respin; /**< if positive, repeat count so total photon=totalphoton*respin; if negative, total number is processed in respin subset */ - unsigned int srcnum; /**< number of sources for simultaneous pattern sources */ - unsigned int savedetflag; /**< number of sources for simultaneous pattern sources */ - int reserved[2]; /**< reserved fields for future extension */ + char magic[4]; // magic bits= 'M','C','X','H' + unsigned int version; // version of the mch file format + unsigned int maxmedia; // number of media in the simulation + unsigned int detnum; // number of detectors in the simulation + unsigned int colcount; // how many output files per detected photon + unsigned int totalphoton; // how many total photon simulated + unsigned int detected; // how many photons are detected (not necessarily all saved) + unsigned int savedphoton; // how many detected photons are saved in this file + float unitinmm; // what is the voxel size of the simulation + unsigned int seedbyte; // how many bytes per RNG seed + float normalizer; // what is the normalization factor + int respin; // if positive, repeat count so total photon=totalphoton*respin; if negative, total number is processed in respin subset + unsigned int srcnum; // number of sources for simultaneous pattern sources + unsigned int savedetflag; // number of sources for simultaneous pattern sources + unsigned int totalsource; // total source number when multiple sources are defined + int reserved[1]; // reserved fields for future extension } History; When the `-q` flag is set to 1, the detected photon initial seeds are also stored @@ -1360,29 +1361,29 @@ a .jdat file is written, which is a pure JSON file. This file contains a hierach record of the following JSON structure { - "MCXData": { + "MCXData":{ "Info":{ "Version": - "MediaNum": - "DetNum": - ... - "Media":{ - ... - } + "MediaNum": + "DetNum": + ... + "Media":{ + ... + } }, "PhotonData":{ "detid": - "nscat": - "ppath": - "mom": - "p": - "v": - "w0": + "nscat": + "ppath": + "mom": + "p": + "v": + "w0": }, "Trajectory":{ "photonid": - "p": - "w0": + "p": + "w0": }, "Seed":[ ... diff --git a/pmcx/Makefile b/pmcx/Makefile new file mode 100644 index 00000000..ba2eef21 --- /dev/null +++ b/pmcx/Makefile @@ -0,0 +1,6 @@ +all: + python3 -m pip wheel . +clean: + rm -rf *.whl +install: + python3 -m pip install --force-reinstall pmcx-*.whl diff --git a/pmcx/README.md b/pmcx/README.md index 5dd0c814..63600c3a 100644 --- a/pmcx/README.md +++ b/pmcx/README.md @@ -4,7 +4,7 @@ - Copyright: (C) Matin Raayai Ardakani (2022-2023) , Qianqian Fang (2019-2023) , Fan-Yu Yen (2023) - License: GNU Public License V3 or later -- Version: 0.2.10 +- Version: 0.2.11 - URL: https://pypi.org/project/pmcx/ - Github: https://github.com/fangq/mcx diff --git a/pmcx/pmcx/__init__.py b/pmcx/pmcx/__init__.py index 21c3844d..7b4ae323 100644 --- a/pmcx/pmcx/__init__.py +++ b/pmcx/pmcx/__init__.py @@ -49,7 +49,7 @@ # from .files import loadmc2, loadmch, load, save from .bench import bench -__version__ = "0.2.10" +__version__ = "0.2.11" __all__ = ( "gpuinfo", diff --git a/pmcx/pmcx/utils.py b/pmcx/pmcx/utils.py index 1c691620..d8e87f1c 100644 --- a/pmcx/pmcx/utils.py +++ b/pmcx/pmcx/utils.py @@ -357,6 +357,7 @@ def detphoton(detp, medianum, savedetflag, issaveref=None, srcnum=None): newdetp['p'] or ['v']: exit position and direction, when cfg.issaveexit=1 (3) newdetp['w0']: photon initial weight at launch time (3) newdetp['s']: exit Stokes parameters for polarized photon (4) + newdetp['srcid']: the ID of the source when multiple sources are defined (1) """ newdetp = {} c0 = 0 @@ -367,6 +368,9 @@ def detphoton(detp, medianum, savedetflag, issaveref=None, srcnum=None): newdetp["w0"] = detp[0, :].transpose() else: newdetp["detid"] = detp[0, :].astype(int).transpose() + if np.any(newdetp["detid"] > 65535): + newdetp["srcid"] = np.right_shift(newdetp["detid"], 16, dtype=np.int32) + newdetp["detid"] = np.bitwise_and(newdetp["detid"], 0xff, dtype=np.int32) c0 = 1 length = medianum diff --git a/pmcx/setup.py b/pmcx/setup.py index fd04c75c..0e2adb7d 100644 --- a/pmcx/setup.py +++ b/pmcx/setup.py @@ -123,7 +123,7 @@ def build_extension(self, ext): setup( name="pmcx", packages=["pmcx"], - version="0.2.10", + version="0.2.11", requires=["numpy"], license="GPLv3+", author="Matin Raayai Ardakani, Qianqian Fang, Fan-Yu Yen", diff --git a/src/mcx_core.cu b/src/mcx_core.cu index 4abbb5e9..6e6bd0c7 100644 --- a/src/mcx_core.cu +++ b/src/mcx_core.cu @@ -374,11 +374,11 @@ __device__ inline void savedetphoton(float n_det[], uint* detectedphoton, float* baseaddr *= gcfg->reclen; if (SAVE_DETID(gcfg->savedetflag)) { - n_det[baseaddr++] = detid; - if (gcfg->extrasrclen && gcfg->srcid <= 0) { - n_det[baseaddr++] = (((int)ppath[gcfg->w0offset - 1]) << 16); + detid |= (((int)ppath[gcfg->w0offset - 1]) << 16); } + + n_det[baseaddr++] = detid; } for (i = 0; i < gcfg->partialdata; i++) { diff --git a/src/mcx_utils.c b/src/mcx_utils.c index 5607e048..b85a1ad9 100644 --- a/src/mcx_utils.c +++ b/src/mcx_utils.c @@ -1024,12 +1024,13 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c } if (cfg->his.detected == 0 && cfg->his.savedphoton) { - char colnum[] = {1, 3, 1}; - char* dtype[] = {"uint32", "single", "single"}; - char* dname[] = {"photonid", "p", "w0"}; + char colnum[] = {1, 3, 1, 1}; + char* dtype[] = {"uint32", "single", "single", "uint32"}; + char* dname[] = {"photonid", "p", "w0", "srcid"}; + int activecol = sizeof(colnum) - (cfg->his.totalsource == 1); cJSON_AddItemToObject(obj, "Trajectory", dat = cJSON_CreateObject()); - for (int id = 0; id < sizeof(colnum); id++) { + for (int id = 0; id < activecol; id++) { uint dims[2] = {count, colnum[id]}; float* buf = (float*)calloc(dims[0] * dims[1], sizeof(float)); @@ -1048,9 +1049,9 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c col += dims[1]; } } else { - char colnum[] = {1, cfg->his.maxmedia, cfg->his.maxmedia, cfg->his.maxmedia, 3, 3, 1}; - char* dtype[] = {"uint32", "uint32", "single", "single", "single", "single", "single"}; - char* dname[] = {"detid", "nscat", "ppath", "mom", "p", "v", "w0"}; + char colnum[] = {1, cfg->his.maxmedia, cfg->his.maxmedia, cfg->his.maxmedia, 3, 3, 1, 4}; + char* dtype[] = {"uint32", "uint32", "single", "single", "single", "single", "single", "single"}; + char* dname[] = {"detid", "nscat", "ppath", "mom", "p", "v", "w0", "s"}; cJSON_AddItemToObject(obj, "PhotonData", dat = cJSON_CreateObject()); for (int id = 0; id < sizeof(colnum); id++) { @@ -1059,6 +1060,7 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c void* val = NULL; float* fbuf = NULL; uint* ibuf = NULL; + int hassrcid = ((cfg->his.totalsource > 1) && id == 1); if (!strcmp(dtype[id], "uint32")) { ibuf = (uint*)calloc(dims[0] * dims[1], sizeof(uint)); @@ -1066,6 +1068,10 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c for (int i = 0; i < dims[0]; i++) for (int j = 0; j < dims[1]; j++) { ibuf[i * dims[1] + j] = ppath[i * cfg->his.colcount + col + j]; + + if (hassrcid) { + ibuf[i * dims[1] + j] &= 0xFFFF; + } } val = (void*)ibuf; @@ -1080,6 +1086,23 @@ void mcx_savejdet(float* ppath, void* seeds, uint count, int doappend, Config* c val = (void*)fbuf; } + if (hassrcid) { + uint* srcid = (uint*)calloc(dims[0] * dims[1], sizeof(uint)); + + for (int i = 0; i < dims[0] * dims[1]; i++) { + srcid[i] = ibuf[i] >> 16; + ibuf[i] &= 0xFFFF; + } + + cJSON_AddItemToObject(dat, "srcid", sub = cJSON_CreateObject()); + + if (mcx_jdataencode((void*)srcid, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, 0, cfg)) { + MCX_ERROR(-1, "error when converting to JSON"); + } + + free(srcid); + } + cJSON_AddItemToObject(dat, dname[id], sub = cJSON_CreateObject()); if (mcx_jdataencode(val, 2, dims, dtype[id], 4, cfg->zipid, sub, 0, 0, cfg)) { diff --git a/utils/loadmch.m b/utils/loadmch.m index 0fdb5d2e..ccf00709 100644 --- a/utils/loadmch.m +++ b/utils/loadmch.m @@ -28,7 +28,7 @@ % [detid(1) nscat(M) ppath(M) mom(M) p(3) v(3) w0(1) s(4)] % header: file header info, a structure has the following fields % [version,medianum,detnum,recordnum,totalphoton,detectedphoton, -% savedphoton,lengthunit,seedbyte,normalizer,respin,srcnum,savedetflag] +% savedphoton,lengthunit,seedbyte,normalizer,respin,srcnum,savedetflag,totalsource] % photonseed: (optional) if the mch file contains a seed section, this % returns the seed data for each detected photon. Each row of % photonseed is a byte array, which can be used to initialize a @@ -70,7 +70,8 @@ respin=fread(fid,1,'int'); srcnum=fread(fid,1,'uint'); savedetflag=fread(fid,1,'uint'); - junk=fread(fid,2,'uint'); + totalsource=fread(fid,1,'uint'); + junk=fread(fid,1,'uint'); detflag=dec2bin(bitand(savedetflag,(2^8-1)))-'0'; if(strcmp(endian,'ieee-le')) @@ -119,5 +120,5 @@ 'recordnum',header(4),'totalphoton',header(5),... 'detectedphoton',header(6),'savedphoton',header(7),... 'lengthunit',header(8),'seedbyte',seedbyte,'normalizer',normalizer,... - 'respin',respin,'srcnum',srcnum,'savedetflag',savedetflag); + 'respin',respin,'srcnum',srcnum,'savedetflag',savedetflag,'totalsource',totalsource); end diff --git a/utils/mcxdetphoton.m b/utils/mcxdetphoton.m index 4d96014c..c904ff1d 100644 --- a/utils/mcxdetphoton.m +++ b/utils/mcxdetphoton.m @@ -38,6 +38,10 @@ newdetp.w0=detp(1,:)'; else newdetp.detid=int32(detp(1,:))'; + if(any(newdetp.detid>65535)) + newdetp.srcid=bitshift(newdetp.detid,-16); + newdetp.detid=bitand(newdetp.detid,int32(hex2dec('ffff'))); + end end c0=2; end