Skip to content

Commit

Permalink
remove redundant functions in mcxlab
Browse files Browse the repository at this point in the history
  • Loading branch information
@fangq
fangq committed Sep 7, 2023
1 parent 53e2681 commit 24c3533
Show file tree
Hide file tree
Showing 6 changed files with 47 additions and 210 deletions.
7 changes: 4 additions & 3 deletions mcxlab/examples/demo_fullhead_atlas.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,13 +5,14 @@
% atlas template(USC 19.5 year group[Sanchez2012]).
%
% This demo is identical to the MCX simulation used for Fig.9(a) in
% TranYan2019(submitted).
% TranYan2020.
%
% [Sanchez2012] C.E.Sanchez J.E.Richards and C.R.Almli, “Age-Specific MRI Templates
% for Pediatric Neuroimaging,” Developmental Neuropsychology 37, 379–399 (2012).
%
% [TranYan2019] A.P.Tran, S.Yan and Q.Fang, "Improving model-based fNIRS
% analysis using mesh-based anatomical and light-transport models".
% [TranYan2020] Tran AP+, Yan S+, Fang Q*, (2020) “Improving model-based fNIRS
% analysis using mesh-based anatomical and light-transport models,"
% Neurophotonics, 7(1), 015008
%
% This file is part of Monte Carlo eXtreme (MCX) URL:http://mcx.sf.net
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Expand Down
2 changes: 1 addition & 1 deletion mcxlab/examples/demo_replay_frequencydomain.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@

% MCXLAB SETTINGS
clear cfg
cfg.nphoton=1e9;
cfg.nphoton=3e7;
cfg.vol=uint8(layered_model);
cfg.prop=opt_params;
cfg.tstart=0;
Expand Down
36 changes: 26 additions & 10 deletions mcxlab/mcxlab.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,29 +43,39 @@
% the 2D plane in such case.
% for 2D simulations, Example: <demo_mcxlab_2d.m>
%
% MCXLAB also accepts 4D arrays to define continuously varying media.
% MCXLAB also accepts 4D arrays to define continuously varying media.
% The following formats are accepted
% 1 x Nx x Ny x Nz float32 array: mua values for each voxel (must use permute to make 1st dimension singleton)
% 2 x Nx x Ny x Nz float32 array: mua/mus values for each voxel (g/n use prop(2,:))
% 4 x Nx x Ny x Nz uint8 array: mua/mus/g/n gray-scale (0-255) interpolating between prop(2,:) and prop(3,:)
% 2 x Nx x Ny x Nz uint16 array: mua/mus gray-scale (0-65535) interpolating between prop(2,:) and prop(3,:)
% 1 x Nx x Ny x Nz float32 array: "mua_float"/101 format: mua values for each voxel (must use permute to make 1st dimension singleton)
% 2 x Nx x Ny x Nz float32 array: "muamus_float"/100 format: mua/mus values for each voxel (g/n use prop(2,:))
% 3 x Nx x Ny x Nz float32 array: "labelplus"/99 format: {[float32: property_value][float32: property_index i, 0-3][float32: tissue label]},
% the optical properties are first read based on the label, then the i-th property is replaced by the property value
% 3 x Nx x Ny x Nz int16/uint16 array: "mixlabel"/98 format: {[int16 label1][uint16 label2][0-65535 label1 percentage]}
% 4 x Nx x Ny x Nz uint8 array: "asgn_byte"/103 format: mua/mus/g/n gray-scale (0-255) interpolating between prop(2,:) and prop(3,:)
% 2 x Nx x Ny x Nz uint16 array: "muamus_short"/104 format: mua/mus gray-scale (0-65535) interpolating between prop(2,:) and prop(3,:)
% 8 x Nx x Ny x Nz uint8 array: "svmc"/97 format: split-voxel MC (SVMC) hybrid domain, can be created by mcxsvmc.m
% Example: <demo_continuous_mua_mus.m>. If voxel-based media are used, partial-path/momentum outputs are disabled
% Example: <demo_svmc_*.m> for SVMC related examples
% *cfg.prop: an N by 4 array, each row specifies [mua, mus, g, n] in order.
% the first row corresponds to medium type 0
% (background) which is typically [0 0 1 1]. The
% second row is type 1, and so on. The background
% medium (type 0) has special meanings: a photon
% terminates when moving from a non-zero to zero voxel.
% medium (type 0) has special meanings: a photon is terminated
% when moving from a non-zero to zero voxel.
% *cfg.tstart: starting time of the simulation (in seconds)
% *cfg.tstep: time-gate width of the simulation (in seconds)
% *cfg.tend: ending time of the simulation (in second)
% *cfg.srcpos: a 1 by 3 vector, the position of the source in grid unit
% *cfg.srcpos: a 1 by 3 vector, the position of the source in grid unit; if a non-zero
% 4th element is given, it specifies the initial weight of each photon packet
% (or a multiplier in the cases of pattern/pattern3d sources); this initial weight
% can be negative - the output fluence is expected to be linearly proportional
% to this initial weight (w0) before normalization.
% *cfg.srcdir: a 1 by 3 vector, specifying the incident vector; if srcdir
% contains a 4th element, it specifies the focal length of
% the source (only valid for focuable src, such as planar, disk,
% fourier, gaussian, pattern, slit, etc); if the focal length
% is nan, all photons will be launched isotropically regardless
% of the srcdir direction.
% of the srcdir direction; if the focal length is -inf, the launch
% angle will be computed based on the Lambertian (cosine) distribution.
%
%== MC simulation settings ==
% cfg.seed: seed for the random number generator (integer) [0]
Expand All @@ -92,7 +102,9 @@
% '0': this face is not used to detector photons
% '1': this face is used to capture photons (if output detphoton)
% see <demo_bc_det.m>
% cfg.isnormalized:[1]-normalize the output fluence to unitary source, 0-no reflection
% cfg.isnormalized:[1]-normalize the output fluence to unitary source, 0-no reflection.
% setting isnormalized to 2 in the replay mode builds the Jacobian
% with Born approximation instead of the default Rytov approximation
% cfg.isspecular: 1-calculate specular reflection if source is outside, [0] no specular reflection
% cfg.maxgate: the num of time-gates per simulation
% cfg.minenergy: terminate photon when weight less than this level (float) [0.0]
Expand Down Expand Up @@ -170,6 +182,10 @@
% set by srcparam1(1) (in grid unit); if srcparam1(2) is set to a non-zero
% value, this source defines a ring (annulus) shaped source, with
% srcparam1(2) denoting the inner circle's radius, here srcparam1(1)>=srcparam1(2)
% 'ring' [*] - a uniform ring/ring-sector source pointing along srcdir; the outer radius
% of the ring is srcparam1(1) (in grid unit) and the inner radius is srcparam1(2);
% srcparam1(3) and srcparam1(4) can optionally set the lower and upper angular
% bound (in positive rad) of a ring sector if at least one of those is positive.
% 'fourierx' [*] - a general Fourier source, the parameters are
% srcparam1: [v1x,v1y,v1z,|v2|], srcparam2: [kx,ky,phi0,M]
% normalized vectors satisfy: srcdir cross v1=v2
Expand Down
18 changes: 12 additions & 6 deletions src/mcx_utils.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
#include <math.h>
#include <ctype.h>
#include <errno.h>
#include <time.h>

#ifndef WIN32
#include <sys/ioctl.h>
Expand Down Expand Up @@ -3344,10 +3345,11 @@ void mcx_replayinit(Config* cfg, float* detps, int dimdetps[2], int seedbyte) {

for (i = 0; i < dimdetps[1]; i++) {
if (cfg->replaydet <= 0 || cfg->replaydet == (int) (detps[i * dimdetps[0]])) {
if (i != cfg->nphoton)
if (i != cfg->nphoton) {
memcpy((char*) (cfg->replay.seed) + cfg->nphoton * seedbyte,
(char*) (cfg->replay.seed) + i * seedbyte,
seedbyte);
}

cfg->replay.weight[cfg->nphoton] = 1.f;
cfg->replay.tof[cfg->nphoton] = 0.f;
Expand Down Expand Up @@ -3444,7 +3446,11 @@ void mcx_validatecfg(Config* cfg, float* detps, int dimdetps[2], int seedbyte) {
cfg->srcpos.z--; /*convert to C index, grid center*/
}

if (cfg->tstart > cfg->tend || cfg->tstep == 0.f) {
if (cfg->tstep == 0.f) {
cfg->tstep = cfg->tend;
}

if (cfg->tstart >= cfg->tend || cfg->tstep == 0.f) {
MCX_ERROR(-6, "incorrect time gate settings");
}

Expand All @@ -3456,10 +3462,6 @@ void mcx_validatecfg(Config* cfg, float* detps, int dimdetps[2], int seedbyte) {
MCX_ERROR(-6, "field 'steps' can not have zero elements");
}

if (cfg->tend <= cfg->tstart) {
MCX_ERROR(-6, "field 'tend' must be greater than field 'tstart'");
}

gates = (int) ((cfg->tend - cfg->tstart) / cfg->tstep + 0.5);

if (cfg->maxgate > gates) {
Expand Down Expand Up @@ -3490,6 +3492,10 @@ void mcx_validatecfg(Config* cfg, float* detps, int dimdetps[2], int seedbyte) {
}
}

if (cfg->seed < 0 && cfg->seed != SEED_FROM_FILE) {
cfg->seed = time(NULL);
}

if ((cfg->outputtype == otJacobian || cfg->outputtype == otWP || cfg->outputtype == otDCS || cfg->outputtype == otRF)
&& cfg->seed != SEED_FROM_FILE) {
MCX_ERROR(-6, "Jacobian output is only valid in the reply mode. Please define cfg.seed");
Expand Down
188 changes: 1 addition & 187 deletions src/mcxlab.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,6 @@
#include <stdio.h>
#include <string.h>
#include <exception>
#include <time.h>
#include <math.h>

#include "mex.h"
Expand Down Expand Up @@ -76,7 +75,6 @@
typedef mwSize dimtype;

void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cfg);
void mcx_validate_config(Config* cfg);
void mcxlab_usage();

float* detps = NULL; //! buffer to receive data from cfg.detphotons field
Expand Down Expand Up @@ -236,7 +234,7 @@ void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) {
cfg.issaveseed = (nlhs >= 4); /** save detected photon seeds to the 4th output if present */

/** Validate all input fields, and warn incompatible inputs */
mcx_validate_config(&cfg);
mcx_validatecfg(&cfg, detps, dimdetps, seedbyte);

partialdata = (cfg.medianum - 1) * (SAVE_NSCAT(cfg.savedetflag) + SAVE_PPATH(cfg.savedetflag) + SAVE_MOM(cfg.savedetflag));
hostdetreclen = partialdata + SAVE_DETID(cfg.savedetflag) + 3 * (SAVE_PEXIT(cfg.savedetflag) + SAVE_VEXIT(cfg.savedetflag)) + SAVE_W0(cfg.savedetflag) + 4 * SAVE_IQUV(cfg.savedetflag);
Expand Down Expand Up @@ -1151,190 +1149,6 @@ void mcx_set_field(const mxArray* root, const mxArray* item, int idx, Config* cf
}
}

/**
* @brief Pre-computing the detected photon weight and time-of-fly from partial path input for replay
*
* When detected photons are replayed, this function recalculates the detected photon
* weight and their time-of-fly for the replay calculations.
*
* @param[in,out] cfg: the simulation configuration structure
*/

void mcx_replay_prep(Config* cfg) {
int i, j, hasdetid = 0, offset;
float plen;

if (cfg->seed == SEED_FROM_FILE && detps == NULL) {
mexErrMsgTxt("you give cfg.seed for replay, but did not specify cfg.detphotons.\nPlease define it as the detphoton output from the baseline simulation");
}

if (detps == NULL || cfg->seed != SEED_FROM_FILE) {
return;
}

if (cfg->nphoton != dimdetps[1]) {
mexErrMsgTxt("the column numbers of detphotons and seed do not match");
}

if (seedbyte == 0) {
mexErrMsgTxt("the seed input is empty");
}

hasdetid = SAVE_DETID(cfg->savedetflag);
offset = SAVE_NSCAT(cfg->savedetflag) * (cfg->medianum - 1);

if (((!hasdetid) && cfg->detnum > 1) || !SAVE_PPATH(cfg->savedetflag)) {
mexErrMsgTxt("please rerun the baseline simulation and save detector ID (D) and partial-path (P) using cfg.savedetflag='dp' ");
}

cfg->replay.weight = (float*)malloc(cfg->nphoton * sizeof(float));
cfg->replay.tof = (float*)calloc(cfg->nphoton, sizeof(float));
cfg->replay.detid = (int*)calloc(cfg->nphoton, sizeof(int));

cfg->nphoton = 0;

for (i = 0; i < dimdetps[1]; i++) {
if (cfg->replaydet <= 0 || cfg->replaydet == (int)(detps[i * dimdetps[0]])) {
if (i != cfg->nphoton) {
memcpy((char*)(cfg->replay.seed) + cfg->nphoton * seedbyte, (char*)(cfg->replay.seed) + i * seedbyte, seedbyte);
}

cfg->replay.weight[cfg->nphoton] = 1.f;
cfg->replay.tof[cfg->nphoton] = 0.f;
cfg->replay.detid[cfg->nphoton] = (hasdetid) ? (int)(detps[i * dimdetps[0]]) : 1;

for (j = hasdetid; j < cfg->medianum - 1 + hasdetid; j++) {
plen = detps[i * dimdetps[0] + offset + j];
cfg->replay.weight[cfg->nphoton] *= expf(-cfg->prop[j - hasdetid + 1].mua * plen);
plen *= cfg->unitinmm;
cfg->replay.tof[cfg->nphoton] += plen * R_C0 * cfg->prop[j - hasdetid + 1].n;
}

if (cfg->replay.tof[cfg->nphoton] < cfg->tstart || cfg->replay.tof[cfg->nphoton] > cfg->tend) { /*need to consider -g*/
continue;
}

cfg->nphoton++;
}
}

cfg->replay.weight = (float*)realloc(cfg->replay.weight, cfg->nphoton * sizeof(float));
cfg->replay.tof = (float*)realloc(cfg->replay.tof, cfg->nphoton * sizeof(float));
cfg->replay.detid = (int*)realloc(cfg->replay.detid, cfg->nphoton * sizeof(int));
}

/**
* @brief Validate all input fields, and warn incompatible inputs
*
* Perform self-checking and raise exceptions or warnings when input error is detected
*
* @param[in,out] cfg: the simulation configuration structure
*/

void mcx_validate_config(Config* cfg) {
int i, gates, idx1d, isbcdet = 0;
const char boundarycond[] = {'_', 'r', 'a', 'm', 'c', '\0'};
const char boundarydetflag[] = {'0', '1', '\0'};
unsigned int partialdata = (cfg->medianum - 1) * (SAVE_NSCAT(cfg->savedetflag) + SAVE_PPATH(cfg->savedetflag) + SAVE_MOM(cfg->savedetflag));
unsigned int hostdetreclen = partialdata + SAVE_DETID(cfg->savedetflag) + 3 * (SAVE_PEXIT(cfg->savedetflag) + SAVE_VEXIT(cfg->savedetflag)) + SAVE_W0(cfg->savedetflag);
hostdetreclen += cfg->polmedianum ? (4 * SAVE_IQUV(cfg->savedetflag)) : 0; // for polarized photon simulation

if (!cfg->issrcfrom0) {
cfg->srcpos.x--;
cfg->srcpos.y--;
cfg->srcpos.z--; /*convert to C index, grid center*/
}

if (cfg->tstart > cfg->tend || cfg->tstep == 0.f) {
mexErrMsgTxt("incorrect time gate settings");
}

if (ABS(cfg->srcdir.x * cfg->srcdir.x + cfg->srcdir.y * cfg->srcdir.y + cfg->srcdir.z * cfg->srcdir.z - 1.f) > 1e-5) {
mexErrMsgTxt("field 'srcdir' must be a unitary vector");
}

if (cfg->steps.x == 0.f || cfg->steps.y == 0.f || cfg->steps.z == 0.f) {
mexErrMsgTxt("field 'steps' can not have zero elements");
}

if (cfg->tend <= cfg->tstart) {
mexErrMsgTxt("field 'tend' must be greater than field 'tstart'");
}

gates = (int)((cfg->tend - cfg->tstart) / cfg->tstep + 0.5);

if (cfg->maxgate > gates) {
cfg->maxgate = gates;
}

if (cfg->sradius > 0.f) {
cfg->crop0.x = MAX((int)(cfg->srcpos.x - cfg->sradius), 0);
cfg->crop0.y = MAX((int)(cfg->srcpos.y - cfg->sradius), 0);
cfg->crop0.z = MAX((int)(cfg->srcpos.z - cfg->sradius), 0);
cfg->crop1.x = MIN((int)(cfg->srcpos.x + cfg->sradius), cfg->dim.x - 1);
cfg->crop1.y = MIN((int)(cfg->srcpos.y + cfg->sradius), cfg->dim.y - 1);
cfg->crop1.z = MIN((int)(cfg->srcpos.z + cfg->sradius), cfg->dim.z - 1);
} else if (cfg->sradius == 0.f) {
memset(&(cfg->crop0), 0, sizeof(uint3));
memset(&(cfg->crop1), 0, sizeof(uint3));
} else {
/*
if -R is followed by a negative radius, mcx uses crop0/crop1 to set the cachebox
*/
if (!cfg->issrcfrom0) {
cfg->crop0.x--;
cfg->crop0.y--;
cfg->crop0.z--; /*convert to C index*/
cfg->crop1.x--;
cfg->crop1.y--;
cfg->crop1.z--;
}
}

if (cfg->seed < 0 && cfg->seed != SEED_FROM_FILE) {
cfg->seed = time(NULL);
}

if ((cfg->outputtype == otJacobian || cfg->outputtype == otWP || cfg->outputtype == otDCS || cfg->outputtype == otRF) && cfg->seed != SEED_FROM_FILE) {
mexErrMsgTxt("Jacobian output is only valid in the reply mode. Please define cfg.seed");
}

for (i = 0; i < cfg->detnum; i++) {
if (!cfg->issrcfrom0) {
cfg->detpos[i].x--;
cfg->detpos[i].y--;
cfg->detpos[i].z--; /*convert to C index*/
}
}

if (cfg->shapedata && strstr(cfg->shapedata, ":") != NULL) {
if (cfg->mediabyte > 4) {
mexErrMsgTxt("rasterization of shapes must be used with label-based mediatype");
}

Grid3D grid = {&(cfg->vol), &(cfg->dim), {1.f, 1.f, 1.f}, 0};

if (cfg->issrcfrom0) {
memset(&(grid.orig.x), 0, sizeof(float3));
}

int status = mcx_parse_shapestring(&grid, cfg->shapedata);

if (status) {
mexErrMsgTxt(mcx_last_shapeerror());
}
}

mcx_preprocess(cfg);

cfg->his.maxmedia = cfg->medianum - 1; /*skip medium 0*/
cfg->his.detnum = cfg->detnum;
cfg->his.srcnum = cfg->srcnum;
cfg->his.colcount = hostdetreclen; /*column count=maxmedia+2*/
cfg->his.savedetflag = cfg->savedetflag;
mcx_replay_prep(cfg);
}

/**
* @brief Error reporting function in the mex function, equivallent to mcx_error in binary mode
*
Expand Down
6 changes: 3 additions & 3 deletions utils/json2mcx.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,11 +87,11 @@
mediaclass='uint8';
if(isfield(json.Domain,'MediaFormat'))
idx=find(ismember({'byte','short','integer','muamus_float',...
'mua_float','muamus_half','asgn_byte','muamus_short'},...
'mua_float','muamus_half','asgn_byte','muamus_short','svmc'},...
lower(json.Domain.MediaFormat)));
if(idx)
typebyte=[1,2,4,8,4,4,4,4];
typenames={'uint8','uint16','uint32','single','single','uint16','uint8','uint16'};
typebyte=[1,2,4,8,4,4,4,4,8];
typenames={'uint8','uint16','uint32','single','single','uint16','uint8','uint16','uint8'};
bytelen=typebyte(idx);
mediaclass=typenames{idx};
else
Expand Down

0 comments on commit 24c3533

Please sign in to comment.