Photon hitting density calculations

[Jesse Lam]

Hello,

I am using VTS_MATLAB_v4.10.0Beta and trying to calculate the photon hitting density given a separated source and detector. The issue I am encountering is that when I use the standard diffusion model and the Monte Carlo method, they appear to give very different results. In particular, for the Monte Carlo approach, the detector (on the right set at 20 mm) seems oddly illuminated. After trying to debug this issue, I would like to look into what "Vts.Factories.ComputationFactory.GetPHD" is doing, but I cannot find where this function is located. Perhaps I am passing the inputs incorrectly?

My parameters are:
Source detector separation: 20 mm

Absorption coefficient: 0.01 mm^-1

Reduced scattering coefficient: 1.6 mm^-1

g: 0.8

n=1.4

# of photons (for the Monte Carlo): 100000

Any help at all would be extremely appreciated.

Thank you!

Kind regards,

Jesse Lam

[Carole Hayakawa]

Hi Jesse,

The PHD is created by using a method of Schotland:

Schotland, J.C., Haselgrove, J.C., and Leigh, J.S., ”Photon Hitting Density” Applied Optics, 32(4):448-453, 1993.

Briefly, it determines the Green's function solution from the source to a point in the tissue and then another Green's from that point to the detector, and combined the two. If you would like to see the code that does this, you can pull a copy of our source code by following instructions here:

https://github.com/VirtualPhotonics/VTS/wiki/Getting-Started-in-Windows

there are instructions for Linux and Mac in case you're not working on Windows.

Once you have a copy of the source code, you can look at a standard diffusion solution, for example,

vts/src/Modeling/ForwardSolvers/PointSourceSDAForwardSolver.cs

and then the method GetPHD to see the code.

I'm a little confused because you say you are working with the MATLAB interop code and that code only runs standard diffusion approximation solutions and *scaled* Monte Carlo solutions (no running of photons here, results are derived by scaling an already run simulation).  So not sure how you are running a *conventional* MC (with photons) and obtaining the PHD using the MATLAB code.

If you'd like please post your MATLAB script and I can take a look.

Best,

Carole

[Jesse Lam]

Thank you for your response and the information, I will see what I can learn from the source code. As for my script, please excuse the mess as it was still a work-in-progress. I have pasted the relevant bits down below. Essentially, I looked at the existing PHD script included with VTS_MATLAB_v4.10.0Beta and replaced the fluence calculation. In other words, instead of the default "fluence=fs.FluenceOfRhoAndZ(...)" I use "fluence=MC_fluence(...)". Within "MC_fluence(...)" I am utilizing "VtsMonteCarlo.RunSimulations(...)". However, as you mentioned, this function is only using a scaled simulation? In that case, perhaps I need to perform the simulation using MCCL and then use that fluence to calculate the PHD?

%% Options
sd=20; %source-detector separation
mua=0.01; %absorption
musp=1.6; %reduced scattering
res=100; %bins
Nphot=100000; %number of photons
g=0.8;
n=1.4;
thick=20;

%% Initialization
% Setting up
op = [mua musp g n];
rhos = linspace(0,sd,res); %detectors
zs = linspace(0,thick,res); %z range

nop = size(op,1);
op_net = NET.createArray('Vts.OpticalProperties', nop);
for i=1:nop
    op_net(i) = Vts.OpticalProperties;
    op_net(i).Mua =  op(i,1);
    op_net(i).Musp = op(i,2);
    op_net(i).G =    op(i,3);
    op_net(i).N =    op(i,4);
end;

fs =  Vts.Modeling.ForwardSolvers.PointSourceSDAForwardSolver();
fluence=MC_fluence(mua,musp,g,n,sd,thick,res,Nphot); %run Monte Carlo

%% Calculating photon hitting density (phd)
phd = Vts.Factories.ComputationFactory.GetPHD(fs, NET.convertArray(fluence,'System.Double'),...
    sd, op_net, NET.convertArray(rhos,'System.Double'), NET.convertArray(zs,'System.Double'));
final = reshape(double(NET.invokeGenericMethod('System.Linq.Enumerable','ToArray',{'System.Double'},phd)),[length(zs) length(rhos) nop]);

Thank you again for all of your helpful advice.

Kind regards,

Jesse

[Carole Hayakawa]

Hi Jesse,

If you'd like to use Monte Carlo to generate PHDs there are methods to do this:

1) Coupled Forward-Adjoint Monte Carlo Simulations of Radiative Transport for the Study of Optical Probe Design in Heterogeneous Tissues.'' SIAM J. on Appl. Math., 68(1), 253-270, 2007.
2) Coupled forward-adjoint Monte Carlo simulation of spatial-angular light fields to determine optical diagnostic sensitivity in turbid media.'' Journal of Biomedical Optics, 19(6), 065003, 2014.

The general idea in these two papers is that the second Green's function from a point in the tissue to the detector (as described above) is replaced by the adjoint Green's function from detector to point in the tissue, and combined with the Green's function from the source to the point in the tissue.  In addition, the angular variation of the light field is used to combine the two Green's functions.

So basically you run a estimate of Radiance as a function of x,y,z and theta,phi from the source.  Then run an estimate of Radiance(x,y,z,theta,phi) from the detector location and match the two spatially and angularly over the voxels within the tissue.  The matching can be done via a MATLAB script.

Let me know if you have more questions.

Carole

[Jesse Lam]

Thank you so much Carole, I'll try my hand at implementing this.