JyangOptics
commented
1 year ago Dear Prof. Fang,
Thanks for your reply. I am trying to use mcxlab to simulate OCT following the method below: [image: image.png]
I've used the mcxdetweight for get the photon weight and used the weighted photon path length as you suggested in the forum: https://groups.google.com/g/mcx-users/c/hZRt9AQx-gc/m/ZXHnaUbx2OMJ
but the results look strange. Could you help me to take a look at the code below:
% OCT simulation of 5 layer structure clear all; close all; % define the stucture geometry cfg.unitinmm = 0.003; % Voxel edge length in mm. use 10 micron here. % define a 5 layer structure cfg.vol=ones(100,100,50); cfg.vol(:,:,5:10)=2; cfg.vol(:,:,11:15)=3; cfg.vol(:,:,16:20)=4; cfg.vol(:,:,21:end)=5; cfg.vol=uint8(cfg.vol); %volumeViewer(cfg.vol); % define the optical properties % format: [mua(1/mm) mus(1/mm) g n] % using the skin data from "Monte-Carlo simulation of OCT structural images of human % skin using experimental B-scans and voxel based approach to optical properties % distribution" cfg.prop=[0 0 1 1 % medium 0: the environment 0.02 35 0.9 1.49 % medium 1: Stratum corneum 0.015 5 0.95 1.37 % medium 2: Upper epidermis 0.02 12 0.85 1.38 % medium 3: Epidermis 0.08 9 0.9 1.36 % medium 4: Dermis 0.2 12 0.95 1.35]; % medium5: Lower Dermis % volumeViewer(cfg.vol); % set seed to make the simulation repeatible cfg.seed=hex2dec('623F9A9E'); % random number generator (RNG) seed. % set number of photon package cfg.nphoton=1e8; % define the source position cfg.srcpos=[50,50,1]; cfg.srcdir=[0 0 1]; % define detector detposition = [50 50 1]; detradius = 1; cfg.detpos = [detposition detradius]; % time-domain simulation parameters cfg.tstart=0; cfg.tend=5e-9; cfg.tstep=5e-10; % GPU thread configuration cfg.autopilot=1; cfg.gpuid=1; %% running simulation with boundary reflection enabled fprintf('running simulation ... this takes about 50 seconds on a GTX 470\n' ); %cfg.isreflect=1; % enable reflection at exterior boundary cfg.isrefint=1; % enable reflection at interior boundary too tic; [f,det]=mcxlab(cfg); toc; cs = squeeze(sum(f.data(:,detposition(2),:,:),4)); cs = cs'; figure, imshow(log10(cs), []); figure, plot(log(cs(:,50))) % weight of detected photons detweight=mcxdetweight(det,det.prop,det.unitinmm); detweight =detweight'; % weighted average of pathlengths for i = 1:length(det.detid) temp1 = 0; temp2 = 0; for j = 1:size(det.ppath,2) temp1 = temp1+det.ppath(i,j)cfg.unitinmmexp(-det.prop(j+1,1)det.ppath(i,j)cfg.unitinmm); temp2 = temp2+exp(-det.prop(j+1,1)det.ppath(i,j)cfg.unitinmm); end optpathlength(i) = temp1/temp2; end % using coherence gating to construct OCT signal lambda = 840e-6; % central wavlength Dlambda = 50e-6; clength = 2log(2)lambda^2/pi/Dlambda; % coherence length z = linspace(0, 0.15, 2048); % reference length I = zeros(1, length(z)); for i = 1:length(z) for j = 1:length(det.detid) %I(i) = I(i)+sqrt(detweight(j))cos(2pi(2z(i)-optpathlength(j))/lambda)exp(-(2z(i)-optpathlength(j))^2/clength^2); I(i) = I(i)+sqrt(detweight(j))exp(-(2z(i)-optpathlength(j))^2/clength^2); end end figure; plot(z, log10(I+1));
Thank you in advance, Jianlong Yang
Qianqian Fang @.***> 于2022年4月26日周二 02:17写道：
you need to call this function
https://github.com/fangq/mcx/blob/master/utils/mcxdetweight.m
see example
https://github.com/fangq/mcx/blob/master/mcxlab/examples/demo_replay_vs_pmc_timedomain.m
— Reply to this email directly, view it on GitHub https://github.com/fangq/mcx/issues/146#issuecomment-1108890574, or unsubscribe https://github.com/notifications/unsubscribe-auth/AQ2JPFLMMQF3UY3Z45NBTD3VG3OUZANCNFSM5UHZQWMA . You are receiving this because you authored the thread.Message ID: @.***>
you need to call this function
https://github.com/fangq/mcx/blob/master/utils/mcxdetweight.m
see example
https://github.com/fangq/mcx/blob/master/mcxlab/examples/demo_replay_vs_pmc_timedomain.m