# VTS_MATLAB_v2.0.0Preview Photon Hitting Density with Two Layers

Sep 29, 2014 at 6:40 PM
One of the updates described in the VTS_MATLAB_v2.0.0Preview documentation is "A 2-layer standard diffusion approximation solution". I was attempting to plot an example of the photon hitting density for a steady-state two layer medium using the demo code in vts_solver_demo.m and am not having success. I tried updating the optical properties to reflect two layers as shown below. I am also not sure how to specify the thickness of the first layer (i believe the second layer is assumed to be infinitely thick). Any help would be welcomed! Thank you.

%% Example PHDOfRhoAndZ
% Evaluate Photon Hitting Density in cylindrical coordinates
% using one set of optical properties and a distributed gaussian source SDA
% solver type.

op = [
[0.01 1 0.8 1.4];
[0.01 10 0.8 1.4];
];
rhos = linspace(0.1,19.9,100); % s-d separation, in mm
zs = linspace(0.1,19.9,100); % z range in mm

VtsSolvers.SetSolverType('DistributedGaussianSourceSDA');
test = VtsSolvers.PHDOfRhoAndZ(op, rhos, zs, 10);

f = figure;
subplot(2,1,1); imagesc(rhos, zs, log(test(:,:,1)));
subplot(2,1,2); imagesc(rhos, zs, log(test(:,:,2)));
Sep 29, 2014 at 10:32 PM
Hi,

Regarding the software to invoke the 2-layer SDA solution. If you look at the bottom of vts_solver_demo.m there are 3 examples of software that call this solution. These examples are for reflectance, but they show the general parameter settings. The SetSolverType would be 'TwoLayerSDA'. I haven't tried generating the PHD with the 2-layer solution using matlab calls yet, but I scanned the code and I think it should work.

Also, at a short course we held this summer, we used the 2-layer solution in our lab exercises. The GUI for the short course can be found at:
http://www.virtualphotonics.org/vts/shortcourse2014/
On the Fluence/Interrogation Panel, you can select the 2-layer and modify the optical properties for each layer. If you change the Layer Height End value for the top layer, be sure to modify the Begin value for the 2nd layer.

Carole
Sep 30, 2014 at 1:05 AM
Oops! I should have looked a little harder. I'll go investigate these examples. Thanks for you help!
Oct 1, 2014 at 12:44 AM
Hi Carole,

I've attempted to make a PHD function, modeled after the one layer function. I updated the solver to "TwoLayerSDAForwardSolver" and pass in the layerThickness. I am getting an error when I attempt to call the fluence method:

the error is:

## No method 'FluenceOfRhoAndZ' with matching signature found for class 'Vts.Modeling.ForwardSolvers.TwoLayerSDAForwardSolver'.

function r = PHDOfRhoAndZTwoLayer(op, rhos, zs, sd,layerThickness)
        %% PHDOfRhoAndZTwoLayer
%   PHDOfRhoAndZTwoLayer(OP, RHOS, ZS, SD, layerThickness)
%
%   OP is an N x 4 matrix of optical properties
%       eg. OP = [[mua1, mus'1, g1, n1]; [mua2, mus'2, g2, n2]; ...];
%   RHO is an 1 x M array of detector locations (in mm)
%       eg. RHO = [1:10];
%   Z is a 1 x M array of z values (in mm)
%       eg. Z = linspace(0.1,19.9,100);
%   SD is the source-detector separation (in mm)
%   layerThickness is thickness of first layer (in mm)

nop = size(op,1);

fs =  Vts.Modeling.ForwardSolvers.TwoLayerSDAForwardSolver;

nLayers = size(op,2);

for i=1:nop
opArray_net = NET.createArray('Vts.IOpticalPropertyRegion', nLayers);
for j=1:nLayers
zRangeNET = Vts.Common.DoubleRange(0,layerThickness);
opNET = Vts.OpticalProperties(op(i,j,1), op(i,j,2), op(i,j,3), op(i,j,4));
tempOpRegion = Vts.Common.LayerOpticalPropertyRegion(zRangeNET, opNET);
opArray_net(j) = tempOpRegion;
end
end;

fluence = double(fs.FluenceOfRhoAndZ(opArray_net,rhos,zs));
%phd = Vts.Factories.ComputationFactory.GetPHD(fs, fluence, sd, op_net, rhos, zs);
phd = Vts.Factories.ComputationFactory.GetPHD(fs, NET.convertArray(fluence,'System.Double'),...
sd, op_net, NET.convertArray(rhos,'System.Double'), NET.convertArray(zs,'System.Double'));

%r = double(NET.invokeGenericMethod('System.Linq.Enumerable','ToArray',{'System.Double'},phd));
r = reshape(double(NET.invokeGenericMethod('System.Linq.Enumerable','ToArray',{'System.Double'},phd)),[length(zs) length(rhos) nop]);
end
--------------------------------------------------------------------------------

Any thoughts/suggestions? Thank you.
Oct 1, 2014 at 8:24 PM
Hi,

We think you are really close to the solution, let me take a look at the solution that was implemented in the GUI and check the types that are being passed to FluenceOfRhoAndZ.

I will get back to you as soon as I have more information.

Thanks,
Lisa
Oct 4, 2014 at 2:00 AM
Hi,

Here is a working solution for the two layer PHD, it has some hard-coded values and expects 2 sets of optical property values. This solution uses the ComputeFluence method and passes the regions (array of IOpticalPropertyRegions) and the solution domain type to calculate the fluence.

I need to check with one of our other developers about passing the layer thickness for the 2nd layer because this should be an infinite value.

Let me know if you have any other questions.

Thanks,
Lisa
        function r = PHDOfRhoAndZTwoLayer(op, rhos, zs, sd, thickness)
%% PHDOfRhoAndZ
%   PHDOfRhoAndZ(OP, RHOS, ZS, SD, THICKNESS)
%
%   OP is an N x 4 matrix of optical properties
%       eg. OP = [[mua1, mus'1, g1, n1]; [mua2, mus'2, g2, n2]; ...];
%   RHO is an 1 x M array of detector locations (in mm)
%       eg. RHO = [1:10];
%   Z is a 1 x M array of z values (in mm)
%       eg. Z = linspace(0.1,19.9,100);
%   SD is the source-detector separation in
%   THICKNESS is the thickness of the tissue

nop = size(op,1);

fs =  Vts.Modeling.ForwardSolvers.TwoLayerSDAForwardSolver();
solutionDomain = Vts.FluenceSolutionDomainType.FluenceOfRhoAndZ;

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;

regions = NET.createArray('Vts.IOpticalPropertyRegion', 2);
regions(1) = Vts.Common.LayerOpticalPropertyRegion(Vts.Common.DoubleRange(0, thickness), op_net(1));
regions(2) = Vts.Common.LayerOpticalPropertyRegion(Vts.Common.DoubleRange(0, thickness), op_net(2));

independentAxes = NET.createArray('Vts.IndependentVariableAxis', 2);
independentAxes(1) = Vts.IndependentVariableAxis.Rho;
independentAxes(2) = Vts.IndependentVariableAxis.Z;
independentValues = NET.createArray('System.Double[]', 2);
independentValues(1) = rhos;
independentValues(2) = zs;
constantValues = NET.createArray('System.Double', 0);

fluence = double(Vts.Factories.ComputationFactory.ComputeFluence(fs, solutionDomain, independentAxes, independentValues, regions, constantValues));
%fluence = double(fs.FluenceOfRhoAndZ(op_net,rhos,zs));

phd = Vts.Factories.ComputationFactory.GetPHD(fs, NET.convertArray(fluence,'System.Double'),...
sd, op_net, NET.convertArray(rhos,'System.Double'), NET.convertArray(zs,'System.Double'));

%r = double(NET.invokeGenericMethod('System.Linq.Enumerable','ToArray',{'System.Double'},phd));
r = reshape(double(NET.invokeGenericMethod('System.Linq.Enumerable','ToArray',{'System.Double'},phd)),[length(zs) length(rhos)]);
end
Oct 4, 2014 at 4:59 PM
Wow! Thanks Lisa! I struggled with that for a while and couldn't get it to work. Nice job!

I just want to mention that the 2-layer SDA fluence solution can produce negative values for certain optical property choices for the two layers. The solution is obtained by using a Hankel transform using a digital filter method for the quadrature points and this can produce the negative values. This means that the PHD may have negative values as well and so cannot be plotted on a log scale unless those negative values are set to something small and >0. We did this in the code that supports the short course GUI (link listed above) so that we could plot on a log scale.

Carole
Oct 7, 2014 at 12:55 AM
Thanks bunches Lisa/Carole - looks like this is functional (hey, I got it to run!) though yes some things appear hard-coded. this is going to be very helpful.

looking forward to the next release!

Albert
Oct 7, 2014 at 1:07 AM
Hi Lisa,

## I think there might be some issues still with the PHDOfRhoAndZTwoLayer function. If I am only looking at a single two-layer system, I believe the call should be:

opsA = [0.01 1 0.8 1.4; 0.01 1 0.8 1.4]; % top/bottom layer OPs case 1
op(1,:,:) = [opsA];

sd = 15; %source detector separation
rhos = linspace(0.1,19.9,100); % s-d separation, in mm
zs = linspace(0.1,19.9,100); % z range in mm

%VtsSolvers.SetSolverType('DistributedGaussianSourceSDA');
layerThickness = 1; % mm

## test = VtsSolvers.PHDOfRhoAndZTwoLayer(op, rhos, zs, sd,layerThickness);

I get an error in the function PHDOfRhoAndZTwoLayer at this line:
regions(2) = Vts.Common.LayerOpticalPropertyRegion(Vts.Common.DoubleRange(0, thickness), op_net(2));

op_net(2) does not exist.

## I think we need some sort of nested loop. The outer loop is over the number of two-layer systems you are running, the inner loop is over the number of layers in your system. This is the convention that the other TwoLayer functions seem to be using. I don't understand the functions well enough to suggest how to actually fix this though.

nLayers = size(op,2);
    for i=1:nop
for j=1:nLayers
opNET = Vts.OpticalProperties(op(i,j,1), op(i,j,2), op(i,j,3), op(i,j,4));
end

## end;

Thank you.
Oct 7, 2014 at 1:12 AM
btw, sorry my messages are sometimes in bold and sometimes tiny. not sure how I am doing this.
Oct 7, 2014 at 2:01 AM
Hi,
I have the following code and it seems to work for me. Please give it a try and let me know if you still have problems.

%% Example PHDOfRhoAndZ for two-layer
% Evaluate Photon Hitting Density in cylindrical coordinates
% using one set of optical properties and a two-layer SDA
% solver type.

ops = [0.01 1 0.8 1.4; 0.01 1 0.8 1.4]; % top/bottom layer OPs
rhos = linspace(0.1,19.9,100); % s-d separation, in mm
zs = linspace(0.1,19.9,100); % z range in mm
layerThickness = 2; % units mm
sourceDetectorSeparation = 10; % units mm

VtsSolvers.SetSolverType('TwoLayerSDA');
test = VtsSolvers.PHDOfRhoAndZTwoLayer(ops, rhos, zs, sourceDetectorSeparation, layerThickness);

f = figure; imagesc(rhos, zs, log10(test));
axis image;
title('Photon Hitting Density of \rho and z');
xlabel('\rho [mm]')
ylabel('z [mm]')
set(f,'Name','PHD of Rho and z');

Carole
Oct 7, 2014 at 7:14 PM
Hello Again,

simple error on my part. some TwoLayer functions like ROfRhoTwoLayer define optical properties as a 3-dimension array

opsA = [0.01 1 0.8 1.4; 0.01 1 0.8 1.4]; % top/bottom layer OPs case 1
opsB = [0.02 1 0.8 1.4; 0.01 1 0.8 1.4]; % top/bottom layer OPs case 2
op(1,:,:) = [opsA];
op(2,:,:) = [opsB];

and this function PHDOfRhoAndZTwoLayer defines it as a 2-dimensional array

ops = [0.01 1 0.8 1.4; 0.01 1 0.8 1.4]; % top/bottom layer OPs

I ran the PHDOfRhoAndZTwoLayer function for several cases to convince myself the two-layer solution was being implemented, and everything seems a-ok. I did however run into some funny business with a few optical property combinations. For case 1, layer 1=layer 2 and the thickness is 1mm. This works ok. For case 2, I change the thickness to 0.9mm and now the layer 1 PHD looks ok, but the layer 2 PHD is NaN. For case 3, I get NaN for the PHD everywhere. Does the two-layer solution breakdown when the photon transport length is greater than the layer thickness? Is the number of photons making it back to the detector just too small?

Thank you :)

Case 1:
op = [ 0.1 1 0.8 1.4; 0.1 1 0.8 1.4]; % top/bottom layer OPs case 1
rhos = linspace(0.1,8,100); % s-d separation, in mm
zs = linspace(0.1,5,100); % z range in mm
layerThickness = 1; % mm
sd = 5;
test = VtsSolvers.PHDOfRhoAndZTwoLayer(op, rhos, zs, sd,layerThickness);

Case 2:
layerThickness = 0.9; % mm

Case 3:
layerThickness = 0.5; % mm
Oct 7, 2014 at 11:52 PM
Hi,
The two layer code applies Kienle's solution: reference http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao-37-4-779
It assumes that the top layer is at least 1 lstar = 1/(mua+mus') thick. The solution is based on an isotropic point source embedded 1 lstar into the tissue WITHIN the top layer. See the reference for more details. So for the optical properties you're trying lstar~=1mm and therefore your top layer thickness needs to be greater than this.
Carole
Oct 8, 2014 at 1:09 AM
I see. Thanks for clearing that up. I'll have a closer look at the Kienle's solution.

Cheers,

Philip