
One of the updates described in the VTS_MATLAB_v2.0.0Preview documentation is "A 2layer standard diffusion approximation solution". I was attempting to plot an example of the photon hitting density for a steadystate 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); % sd 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)));



Developer
Sep 29, 2014 at 9:32 PM

Hi,
Regarding the software to invoke the 2layer 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 2layer 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 2layer 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 2layer 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




Oops! I should have looked a little harder. I'll go investigate these examples. Thanks for you help!




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:
fluence = double(fs.FluenceOfRhoAndZ(opArray_net,rhos,zs));
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 sourcedetector 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.



Coordinator
Oct 1, 2014 at 7: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



Coordinator
Oct 4, 2014 at 1:00 AM

Hi,
Here is a working solution for the two layer PHD, it has some hardcoded 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 sourcedetector 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



Developer
Oct 4, 2014 at 3: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 2layer 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




Thanks bunches Lisa/Carole  looks like this is functional (hey, I got it to run!) though yes some things appear hardcoded. this is going to be very helpful.
looking forward to the next release!
Albert




Hi Lisa,
I think there might be some issues still with the PHDOfRhoAndZTwoLayer function. If I am only looking at a single twolayer 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); % sd 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 twolayer 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.




btw, sorry my messages are sometimes in bold and sometimes tiny. not sure how I am doing this.



Developer
Oct 7, 2014 at 1: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 twolayer
% Evaluate Photon Hitting Density in cylindrical coordinates
% using one set of optical properties and a twolayer 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); % sd 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




Hello Again,
simple error on my part. some TwoLayer functions like ROfRhoTwoLayer define optical properties as a 3dimension 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 2dimensional 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 twolayer solution was being implemented, and everything seems aok. 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 twolayer 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); % sd 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



Developer
Oct 7, 2014 at 10:52 PM

Hi,
The two layer code applies Kienle's solution: reference
http://www.opticsinfobase.org/ao/abstract.cfm?uri=ao374779
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




I see. Thanks for clearing that up. I'll have a closer look at the Kienle's solution.
Cheers,
Philip

