MCX for SFDI simulation

[Yuan Gao]

Why are all the photon paths I detected 0? It feels quite strange. I also set them as tangential in the routine, which feels a bit odd. Could you please take a look for me

Moreover, I set up a light source that should be the same size as the organization, but in the end, it seemed that the edge was gone

Here is my code

clear cfg;

clear all;

%多层组织结构定义

% uint8/uint16/uint32/single/double，每个体素值对应 cfg.prop 的行索引（如 1 对应 prop(2,:)）；

cfg.vol = ones(30, 30, 50);

cfg.vol(:, :, 20:24) = 2;

cfg.vol(:, :, 25:34) = 3;

cfg.vol(:, :, 35:end) = 4;

cfg.vol = uint8(cfg.vol);

%各层组织光学特性定义

% 介质光学参数，每行对应一种介质：[mua, mus, g, n]：

% - mua：吸收系数（mm⁻¹）；

% - mus：散射系数（mm⁻¹）；

% - g：各向异性因子（-1~1，0 为各向同性散射，0.9 为前向散射）；

% - n：折射率；

% ⚠️ 第 1 行必须是背景介质（类型 0），通常设为 [0, 0, 1, 1]（无吸收、无散射、折射率 1）

cfg.prop = [0 0 1 1 % medium 0: the environment

0.019 7.8 0.89 1.37 % medium 1: skin & skull

0.0004 0.009 0.89 1.37 % medium 2: CSF

0.02 9.0 0.89 1.37 % medium 3: gray matter

0.08 40.9 0.84 1.37]; % medium 4: white matter

% 光源定义

cfg.srctype = 'fourier'; % define an SFDI source

cfg.srcpos = [0 0 49]; % one corner of the illumination area

xphase = pi / 3; % phase offset in the x-dir, must < 2pi

yphase = 0; % phase offset in the x-dir, must < 2pi

%这里的实际fx是kx/L，也就是kx=fx*L

Lx = 30;

Ly = 30;

fx = 0;

fy = 0;

kx = fx*Lx; % wave number in the x-dir

ky = fy*Ly; % wave number in the x-dir

cfg.srcparam1 = [Lx 0 0 kx + xphase / (2 * pi)]; % kx is k-number in x direction

cfg.srcparam2 = [0 Ly 0 ky + yphase / (2 * pi)];

cfg.srcdir = [0 0 -1];

% other simulation information

cfg.nphoton = 3e6;

cfg.seed = 1648335518;

% time-domain simulation parameters

cfg.tstart = 0;

cfg.tend = 5e-9;

cfg.tstep = 5e-10;

% GPU thread configuration

cfg.autopilot = 1;

cfg.gpuid = 1;

% cfg.debuglevel = 'TM';

cfg.issrcfrom0 = 1;

cfg.issaveref=1;

cfg.isreflect = 1;

cfg.issaveexit = 1;

cfg.issavedet = 1;

% cfg.detpos = [15 15 50 2]; % detector position

cfg.detpos=[15 10 49 1;15 20 49 1;10 15 49 1;20 15 49 1];

cfg.vol(:,:,50) = 0;

% calculate the fluence distribution with the given config

[fluence,detpt,vol,seeds,traj]=mcxlab(cfg);

% integrate time-axis (4th dimension) to get CW solutions

cwfluence=sum(fluence.data,4); % fluence rate

cwdref=sum(fluence.dref,4); % diffuse reflectance

% plot configuration and results

subplot(231);

mcxpreview(cfg);title('domain preview');

subplot(232);

imagesc(squeeze(log(cwfluence(:,15,:))));title('fluence at y=15');

subplot(233);

hist(detpt.ppath(:,1),50); title('partial path tissue#1');

subplot(234);

plot(squeeze(fluence.data(15,15,15,:)),'-o');title('TPSF at [15,15,15]');

subplot(235);

newtraj=mcxplotphotons(traj);title('photon trajectories')

subplot(236);

imagesc(squeeze(log(cwdref(:,:,50))));title('diffuse refle. at z=0');

[Qianqian Fang]

hi Yuan,

the issue is that you defined your layered brain model with the +z axis pointing to the scalp-to-white-matter direction, but you placed your source and detectors on the maximum z-plane, making them direct in contact with the white matter instead of the scalp.

to fix it, you need to move them to the bottom of the domain, if you add

cfg.srcpos = [0 0 1];
cfg.srcdir=[0 0 1];
cfg.detpos(:,3)=1;
cfg.vol(:,:,1)=0;  % in order to record diffuse reflectance at z=1 plane

and rerun your simulation, you should see expected results (attached)

for the final subplot, you need to change your cwdref's z index from 50 to 1

imagesc(squeeze(log(cwdref(:,:,1))));title('diffuse refle. at z=0');

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/fa42add7-8567-4632-9ab7-c7efe01d212an%40googlegroups.com.