Subpixel simulation

[Vijitha Periyasamy]

Dear Dr Fang,

I am facing an issue with the subpixel simulation.

The transport mean free path matches the decay constant estimated from the mcx output only when the voxel size is 1mm. 

I did look at your scaling in the source code (https://github.com/fangq/mcx/blob/v2025/src/mcx_core.cu#L3736-L3807). The scaling by a constant doesnt change the decay constant.

and the demo_qtest_subpixel.m

I believe the decay constant estimated by fitting an exponential for the energy deposited summed along z should be the same immaterial of the voxel size. Attached is the curves for 0.5 mm voxel and 1 mm voxel.

Here is the code of demo edited to calculate the tmpf and the decay constant

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% MCXLAB - Monte Carlo eXtreme for MATLAB/Octave by Qianqina Fang

%

% In this example, we test the sub-millimeter voxel feature using

% the settings in examples/quicktest (i.e. comparing run_qtest.sh and

% run_grid2x.sh)

%

% This file is part of Monte Carlo eXtreme (MCX) URL:https://mcx.space

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% only clear cfg to avoid accidentally clearing other useful data

clear cfg;

cfg.nphoton = 1e7;

cfg.vol = uint8(ones(60, 60, 60));

cfg.srcpos = [30 30 1];

cfg.srcdir = [0 0 1];

cfg.gpuid = 1;

cfg.autopilot = 1;

cfg.prop = [0 0 1 1; 0.003 7.7 0 1.0];

cfg.tstart = 0;

cfg.tend = 5e-9;

cfg.tstep = 5e-9;

cfg.outputtype = 'energy';

% calculate the flux distribution with the given config

[f1, det1] = mcxlab(cfg);

cfg.vol = uint8(ones(120, 120, 120));

cfg.srcpos = [60 60 1];

cfg.unitinmm = 0.5;

cfg.nphoton = 8e7; % you need to simulate 8x photons to get the same noise

[f2, det2] = mcxlab(cfg);

figure;

subplot(121);

imagesc(squeeze(log(f1(1).data(:, 30, :, 1))));

colorbar;

axis equal;

axis off;

cl = get(gca, 'clim');

subplot(122);

imagesc(squeeze(log(f2(1).data(:, 60, :, 1))));

colorbar;

axis equal;

axis off;

set(gca, 'clim', cl);

tt=squeeze(sum(sum(f1(1).data,2),1));

uu = squeeze(sum(sum(f2(1).data,2),1));

figure,

plot([1:60],(tt));

hold on

plot([1:120]*0.5*0.5,(uu));

[a1] = fit([10:60]'*1, tt(10:60), 'exp1');

[a2] = fit([10:120]'*0.5, uu(10:120), 'exp1');