Skip to content

Commit

Permalink
support trajectory-only mode with debuglevel=T
Browse files Browse the repository at this point in the history
  • Loading branch information
@fangq
fangq committed May 21, 2023
1 parent 0100212 commit 540931d
Show file tree
Hide file tree
Showing 7 changed files with 65 additions and 46 deletions.
44 changes: 26 additions & 18 deletions mcxlab/mcxlab.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -144,7 +144,7 @@
% 'isotropic' - isotropic source, no param needed
% 'cone' - uniform cone beam, srcparam1(1) is the half-angle in radian
% 'gaussian' [*] - a collimated gaussian beam, srcparam1(1) specifies the waist radius (in voxels)
% 'hyperboloid' - a one-sheeted hyperboloid gaussian beam, srcparam1(1) specifies the waist
% 'hyperboloid' [*] - a one-sheeted hyperboloid gaussian beam, srcparam1(1) specifies the waist
% radius (in voxels), srcparam1(2) specifies distance between launch plane and focus,
% srcparam1(3) specifies rayleigh range
% 'planar' [*] - a 3D quadrilateral uniform planar source, with three corners specified
Expand Down Expand Up @@ -258,6 +258,7 @@
% 'R': debug RNG, output fluence.data is filled with 0-1 random numbers
% 'M': return photon trajectory data as the 5th output
% 'P': show progress bar
% 'T': save photon trajectory data only, as the 1st output, disable flux/detp/seeds outputs
% cfg.maxjumpdebug: [10000000|int] when trajectory is requested in the output,
% use this parameter to set the maximum position stored. By default,
% only the first 1e6 positions are stored.
Expand All @@ -271,6 +272,9 @@
% dimensions specified by [size(vol) total-time-gates].
% The content of the array is the normalized fluence at
% each voxel of each time-gate.
%
% when cfg.debuglevel contains 'T', fluence(i).data stores trajectory
% output, see below
% fluence(i).dref is a 4D array with the same dimension as fluence(i).data
% if cfg.issaveref is set to 1, containing only non-zero values in the
% layer of voxels immediately next to the non-zero voxels in cfg.vol,
Expand Down Expand Up @@ -511,22 +515,26 @@
if(exist('newdetpstruct','var'))
varargout{2}=newdetpstruct;
end
end

if(nargout>=5)
for i=1:length(varargout{5})
data=varargout{5}.data;
if(isempty(data))
continue;
end
traj.pos=data(2:4,:).';
traj.id=typecast(data(1,:),'uint32').';
[traj.id,idx]=sort(traj.id);
traj.pos=traj.pos(idx,:);
traj.data=[single(traj.id)' ; data(2:end,idx)];
newtraj(i)=traj;
end
if(exist('newtraj','var'))
varargout{5}=newtraj;
end
if(nargout>=5 || (~isempty(cfg) && isstruct(cfg) && isfield(cfg, 'debuglevel') && ~isempty(regexp(cfg(1).debuglevel, '[tT]', 'once'))))
outputid=5;
if((isfield(cfg, 'debuglevel') && ~isempty(regexp(cfg(1).debuglevel, '[tT]', 'once'))))
outputid=1;
end
end
for i=1:length(varargout{outputid})
data=varargout{outputid}.data;
if(isempty(data))
continue;
end
traj.pos=data(2:4,:).';
traj.id=typecast(data(1,:),'uint32').';
[traj.id,idx]=sort(traj.id);
traj.pos=traj.pos(idx,:);
traj.data=[single(traj.id)' ; data(2:end,idx)];
newtraj(i)=traj;
end
if(exist('newtraj','var'))
varargout{outputid}=newtraj;
end
end
1 change: 1 addition & 0 deletions src/mcx_const.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,6 +60,7 @@
#define MCX_DEBUG_RNG 1 /**< debug flags: 1 - run RNG testing kernel and return RNG numbers */
#define MCX_DEBUG_MOVE 2 /**< debug flags: 2 - save and output photon trajectory data */
#define MCX_DEBUG_PROGRESS 4 /**< debug flags: 4 - print progress bar */
#define MCX_DEBUG_MOVE_ONLY 8 /**< debug flags: 8 - only save photon trajectory data, disable volume and detphoton output */

#define MEDIA_2LABEL_SPLIT 97 /**< media Format: 64bit:{[byte: lower label][byte: upper label][byte*3: reference point][byte*3: normal vector]} */
#define MEDIA_2LABEL_MIX 98 /**< media format: {[int: label1][int: label2][float32: label1 %]} -> 32bit:{[half: label1 %],[byte: label2],[byte: label1]} */
Expand Down
20 changes: 10 additions & 10 deletions src/mcx_core.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1015,7 +1015,7 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*
if (p->w >= 0.f) {
ppath[gcfg->partialdata] += p->w; //< sum all the remaining energy

if (gcfg->debuglevel & MCX_DEBUG_MOVE) {
if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
savedebugdata(p, ((uint)f->ndone) + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata);
}

Expand Down Expand Up @@ -1454,7 +1454,7 @@ __device__ inline int launchnewphoton(MCXpos* p, MCXdir* v, Stokes* s, MCXtime*
f->ndone++;
updateproperty<islabel, issvmc>(prop, *mediaid, t, *idx1d, media, (float3*)p, nuvox, flipdir);

if (gcfg->debuglevel & MCX_DEBUG_MOVE) {
if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
savedebugdata(p, (uint)f->ndone + threadid * gcfg->threadphoton + umin(threadid, gcfg->oddphotons), gdebugdata);
}

Expand Down Expand Up @@ -1801,7 +1801,7 @@ __global__ void mcx_main_loop(uint media[], OutputType field[], float genergy[],
#endif
}

if (gcfg->debuglevel & MCX_DEBUG_MOVE) {
if (gcfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
savedebugdata(&p, (uint)f.ndone + idx * gcfg->threadphoton + umin(idx, gcfg->oddphotons), gdebugdata);
}
}
Expand Down Expand Up @@ -2843,7 +2843,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) {
CUDA_ASSERT(cudaHostGetDevicePointer((int**)&gprogress, (int*)progress, 0));
*progress = 0;

if (cfg->debuglevel & MCX_DEBUG_MOVE) {
if (cfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
CUDA_ASSERT(cudaMalloc((void**) &gdebugdata, sizeof(float) * (debuglen * cfg->maxjumpdebug)));
}

Expand Down Expand Up @@ -3055,7 +3055,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) {

CUDA_ASSERT(cudaMemset(gdetected, 0, sizeof(float)));

if (cfg->debuglevel & MCX_DEBUG_MOVE) {
if (cfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
uint jumpcount = 0;
CUDA_ASSERT(cudaMemcpyToSymbol(gjumpdebug, &jumpcount, sizeof(uint), 0, cudaMemcpyHostToDevice));
}
Expand Down Expand Up @@ -3279,7 +3279,7 @@ void mcx_run_simulation(Config* cfg, GPUInfo* gpu) {
/**
* If '-D M' is specified, we retrieve photon trajectory data and store those to \c cfg.exportdebugdata and \c cfg.debugdatalen
*/
if (cfg->debuglevel & MCX_DEBUG_MOVE) {
if (cfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
uint debugrec = 0;
CUDA_ASSERT(cudaMemcpyFromSymbol(&debugrec, gjumpdebug, sizeof(uint), 0, cudaMemcpyDeviceToHost));
#pragma omp critical
Expand Down Expand Up @@ -3458,7 +3458,7 @@ is more than what your have specified (%d), please use the -H option to specify
*/
#pragma omp master
{
if (cfg->srctype == MCX_SRC_PATTERN && cfg->srcnum > 1) { // post-processing only for multi-srcpattern
if (cfg->issave2pt && cfg->srctype == MCX_SRC_PATTERN && cfg->srcnum > 1) { // post-processing only for multi-srcpattern
srcpw = (float*)calloc(cfg->srcnum, sizeof(float));
energytot = (float*)calloc(cfg->srcnum, sizeof(float));
energyabs = (float*)calloc(cfg->srcnum, sizeof(float));
Expand Down Expand Up @@ -3499,7 +3499,7 @@ is more than what your have specified (%d), please use the -H option to specify
* in joule when cfg.outputtype='fluence', or energy-loss multiplied by mua (1/mm) per voxel
* (joule/mm) when cfg.outputtype='flux' (default).
*/
if (cfg->isnormalized) {
if (cfg->issave2pt && cfg->isnormalized) {
float* scale = (float*)calloc(cfg->srcnum, sizeof(float));
scale[0] = 1.f;
int isnormalized = 0;
Expand Down Expand Up @@ -3634,7 +3634,7 @@ is more than what your have specified (%d), please use the -H option to specify
*/
#ifndef MCX_CONTAINER

if ((cfg->debuglevel & MCX_DEBUG_MOVE) && cfg->parentid == mpStandalone && cfg->exportdebugdata) {
if ((cfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) && cfg->parentid == mpStandalone && cfg->exportdebugdata) {
cfg->his.colcount = MCX_DEBUG_REC_LEN;
cfg->his.savedphoton = cfg->debugdatalen;
cfg->his.totalphoton = cfg->nphoton;
Expand Down Expand Up @@ -3717,7 +3717,7 @@ is more than what your have specified (%d), please use the -H option to specify
CUDA_ASSERT(cudaFree(gsmatrix));
}

if (cfg->debuglevel & MCX_DEBUG_MOVE) {
if (cfg->debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
CUDA_ASSERT(cudaFree(gdebugdata));
}

Expand Down
22 changes: 14 additions & 8 deletions src/mcx_utils.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,7 +142,7 @@ const char outputtype[] = {'x', 'f', 'e', 'j', 'p', 'm', 'r', 'l', '\0'};
* P: show progress bar
*/

const char debugflag[] = {'R', 'M', 'P', '\0'};
const char debugflag[] = {'R', 'M', 'P', 'T', '\0'};

/**
* Recorded fields for detected photons
Expand Down Expand Up @@ -1458,6 +1458,11 @@ void mcx_preprocess(Config* cfg) {
cfg->srcdir.y *= tmp;
cfg->srcdir.z *= tmp;

if (cfg->debuglevel & MCX_DEBUG_MOVE_ONLY) {
cfg->issave2pt = 0;
cfg->issavedet = 0;
}

for (int i = 0; i < 6; i++)
if (cfg->bc[i] && mcx_lookupindex(cfg->bc + i, boundarycond)) {
MCX_ERROR(-4, "unknown boundary condition specifier");
Expand Down Expand Up @@ -5119,17 +5124,17 @@ where possible parameters include (the first value in [*|*] is the default)\n\
-S [1|0] (--save2pt) 1 to save the flux field; 0 do not save\n\
-F [mc2|...] (--outputformat) fluence data output format:\n\
mc2 - MCX mc2 format (binary 32bit float)\n\
jnii - JNIfTI format (http://openjdata.org)\n\
bnii - Binary JNIfTI (http://openjdata.org)\n\
jnii - JNIfTI format (https://neurojson.org)\n\
bnii - Binary JNIfTI (https://neurojson.org)\n\
nii - NIfTI format\n\
hdr - Analyze 7.5 hdr/img format\n\
tx3 - GL texture data for rendering (GL_RGBA32F)\n\
the bnii/jnii formats support compression (-Z) and generate small files\n\
load jnii (JSON) and bnii (UBJSON) files using below lightweight libs:\n\
MATLAB/Octave: JNIfTI toolbox https://github.com/fangq/jnifti, \n\
MATLAB/Octave: JSONLab toolbox https://github.com/fangq/jsonlab, \n\
MATLAB/Octave: JNIfTI toolbox https://github.com/NeuroJSON/jnifti,\n\
MATLAB/Octave: JSONLab toolbox https://github.com/NeuroJSON/jsonlab,\n\
Python: PyJData: https://pypi.org/project/jdata\n\
JavaScript: JSData: https://github.com/fangq/jsdata\n\
JavaScript: JSData: https://github.com/NeuroJSON/jsdata\n\
-Z [zlib|...] (--zip) set compression method if -F jnii or --dumpjson\n\
is used (when saving data to JSON/JNIfTI format)\n\
0 zlib: zip format (moderate compression,fast) \n\
Expand All @@ -5140,7 +5145,7 @@ where possible parameters include (the first value in [*|*] is the default)\n\
5 lz4: LZ4 format (low compression,extrem. fast)\n\
6 lz4hc: LZ4HC format (moderate compression,fast)\n\
--dumpjson [-,0,1,'file.json'] export all settings, including volume data using\n\
JSON/JData (http://openjdata.org) format for \n\
JSON/JData (https://neurojson.org) format for\n\
easy sharing; can be reused using -f\n\
if followed by nothing or '-', mcx will print\n\
the JSON to the console; write to a file if file\n\
Expand All @@ -5157,9 +5162,10 @@ where possible parameters include (the first value in [*|*] is the default)\n\
== Debug options ==\n" S_RESET"\
-D [0|int] (--debug) print debug information (you can use an integer\n\
or or a string by combining the following flags)\n\
-D [''|RMP] 1 R debug RNG\n\
-D [''|RMPT] 1 R debug RNG\n\
/case insensitive/ 2 M store photon trajectory info\n\
4 P print progress bar\n\
8 T save trajectory data only, disable flux/detp\n\
combine multiple items by using a string, or add selected numbers together\n\
\n"S_BOLD S_CYAN"\
== Additional options ==\n" S_RESET"\
Expand Down
15 changes: 8 additions & 7 deletions src/mcxlab.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,7 +177,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) {
if (nlhs >= 1 || (cfg.debuglevel & MCX_DEBUG_MOVE_ONLY)) {
plhs[0] = mxCreateStructMatrix(ncfg, 1, 4, datastruct);
}

Expand Down Expand Up @@ -269,7 +269,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
cfg.seeddata = malloc(cfg.maxdetphoton * sizeof(float) * RAND_WORD_LEN);
}

if (nlhs >= 5) {
if (nlhs >= 5 || (cfg.debuglevel & MCX_DEBUG_MOVE_ONLY)) {
cfg.exportdebugdata = (float*)malloc(cfg.maxjumpdebug * sizeof(float) * MCX_DEBUG_REC_LEN);
cfg.debuglevel |= MCX_DEBUG_MOVE;
}
Expand Down Expand Up @@ -311,15 +311,16 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
fielddim[5] = 1;

/** if 5th output presents, output the photon trajectory data */
if (nlhs >= 5) {
if (nlhs >= 5 || (cfg.debuglevel & MCX_DEBUG_MOVE_ONLY)) {
int outputidx=(cfg.debuglevel & MCX_DEBUG_MOVE_ONLY) ? 0 : 4;
fielddim[0] = MCX_DEBUG_REC_LEN;
fielddim[1] = cfg.debugdatalen; // his.savedphoton is for one repetition, should correct
fielddim[2] = 0;
fielddim[3] = 0;
mxSetFieldByNumber(plhs[4], jstruct, 0, mxCreateNumericArray(2, fielddim, mxSINGLE_CLASS, mxREAL));
mxSetFieldByNumber(plhs[outputidx], jstruct, 0, mxCreateNumericArray(2, fielddim, mxSINGLE_CLASS, mxREAL));

if (cfg.debuglevel & MCX_DEBUG_MOVE) {
memcpy((float*)mxGetPr(mxGetFieldByNumber(plhs[4], jstruct, 0)), cfg.exportdebugdata, fielddim[0]*fielddim[1]*sizeof(float));
if (cfg.debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
memcpy((float*)mxGetPr(mxGetFieldByNumber(plhs[outputidx], jstruct, 0)), cfg.exportdebugdata, fielddim[0]*fielddim[1]*sizeof(float));
}

if (cfg.exportdebugdata) {
Expand Down Expand Up @@ -930,7 +931,7 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
printf("mcx.outputtype='%s';\n", outputstr);
} else if (strcmp(name, "debuglevel") == 0) {
int len = mxGetNumberOfElements(item);
const char debugflag[] = {'R', 'M', 'P', '\0'};
const char debugflag[] = {'R', 'M', 'P', 'T', '\0'};
char debuglevel[MAX_SESSION_LENGTH] = {'\0'};

if (!mxIsChar(item) || len == 0) {
Expand Down
6 changes: 3 additions & 3 deletions src/pmcx.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -910,7 +910,7 @@ py::dict pmcx_interface(const py::dict& user_cfg) {
mcx_config.seeddata = malloc(mcx_config.maxdetphoton * sizeof(float) * RAND_WORD_LEN);
}

if (mcx_config.debuglevel & MCX_DEBUG_MOVE) {
if (mcx_config.debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
mcx_config.exportdebugdata = (float*) malloc(mcx_config.maxjumpdebug * sizeof(float) * MCX_DEBUG_REC_LEN);
}

Expand Down Expand Up @@ -947,14 +947,14 @@ py::dict pmcx_interface(const py::dict& user_cfg) {
field_dim[4] = 1;
field_dim[5] = 1;

if (mcx_config.debuglevel & MCX_DEBUG_MOVE) {
if (mcx_config.debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
field_dim[0] = MCX_DEBUG_REC_LEN;
field_dim[1] = mcx_config.debugdatalen; // his.savedphoton is for one repetition, should correct
field_dim[2] = 0;
field_dim[3] = 0;
auto photon_traj_data = py::array_t<float, py::array::f_style>({field_dim[0], field_dim[1]});

if (mcx_config.debuglevel & MCX_DEBUG_MOVE) {
if (mcx_config.debuglevel & (MCX_DEBUG_MOVE | MCX_DEBUG_MOVE_ONLY)) {
memcpy(photon_traj_data.mutable_data(), mcx_config.exportdebugdata, field_dim[0] * field_dim[1] * sizeof(float));
}

Expand Down
3 changes: 3 additions & 0 deletions utils/mcxplotphotons.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,9 @@
% License: GPLv3, see http://mcx.sf.net for details
%

if(isstruct(traj) && ~isfield(traj, 'id') && isfield(traj, 'data'))
traj=struct('id',typecast(single(traj.data(1,:)'),'uint32'),'pos',traj.data(2:4,:)','weight',traj.data(5,:)', 'data', traj.data);
end
if(~isstruct(traj) && size(traj,2)==6)
traj=struct('id',typecast(single(traj(:,1)),'uint32'),'pos',traj(:,2:4),'weight',traj(:,5));
end
Expand Down

0 comments on commit 540931d

Please sign in to comment.