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');

[Qianqian Fang]

I believe this result is expected.

your script sets the output type to 'energy'. this accumulates the energy-loss within each voxel - if you have smaller voxels, your per-voxel accumulated energy loss will be correspondingly smaller as expected.

your plot is also problematic

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));

while you summed the x/y dimensions, your z-slice thickness also has a 2x difference, which is not accounted for in your plot.

if you set the outputtype to fluence, and plot the fluence over distance, you will see that the two solutions match

cfg.outputtype = 'fluence';
...
tt=squeeze(f1(1).data(30,30,:));
uu =squeeze(f2(1).data(60,60,:));
figure,
semilogy([1:60],(tt));
hold on
semilogy([1:120]*0.5,(uu));

the difference is that each value in the 'energy' output is an integral (over voxels), and each value in fluence/flux output is a point-sample of an underlying distribution - after normalization, the latter is voxel-size independent, but the first one is not.

Qianqian

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/mcx-users/a8647e6c-b47f-450f-9590-e3485899bb85n%40googlegroups.com.